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Abstract 



This thesis presents the original results I have obtained during the three- 
year doctoral period in the Division de Matematicas Aplicadas y Sistemas 
Computacionales (DMASC) of the Instituto Potosino de Investigation Cientffica 
y Tecnologica (IPICYT), in San Luis Potosf, Mexico. These results have been 
obtained under supervision and collaboration of Dr. Haret C. Rosu in what refers 
to the first part of the thesis, and of Dr. Ricardo Femat for the second part. 

The first part deals with some types of factorization methods that we were able to 
develop and that lead us to particular solutions of travelling kink type for reaction- 
diffusion equations and also to more general nonlinear differential equations 
of interest in biology and nonlinear physics. We also applied supersymmetric 
approaches in the context of biological dynamics of microtubules and the related 
transport properties associated to their domain walls. In addition, a complex 
supersymmetric extension of the classical harmonic oscillator by which we obtain 
new oscillatory modes has been developed; results that could be extended to 
physical optics and the physics of cavities. Moreover, an application to chemical 
physics of diatomic molecules using supersymmetric and factorization procedures 
is developed. 

The second part contains a detailed study on the synchronization of the chaotic 
dynamics of two Hodgkin-Huxley neurons, by means of the mathematical tools 
belonging to the geometrical control theory. Despite using different parameters 
for each of the two neurons our analysis shows that synchronization states are 
achieved. The synchronization is attained by the feedback structure of the 
interconnection (coupling). Numerical results for the obtained neuronal dynamical 
states are displayed. 
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Resumen 



Esta tesis presenta los resultados originales que he obtenido durante los tres 
afios de periodo doctoral en la Division de Matematicas Aplicadas y Sistemas 
Computacionales (DMASC) del Instituto Potosino de Investigation Cienti'fica y 
Tecnologica (IPICYT), en San Luis Potosf, Mexico. Estos resultados se han 
obtenido bajo la supervision y colaboracion del Dr. Haret C. Rosu en lo referente 
a la primera parte de tesis, y del Dr. Ricardo Femat para la segunda parte. 

La primera parte trata con algunos metodos de factorizacion que fuimos capaces 
de desarrollar y que nos condujeron a soluciones particulares del tipo kink viajeras 
para ecuaciones de reaccion-difusion y tambien para ecuaciones diferenciales no 
lineales mas generales de interes en biologfa y fisica no lineal. Se aplicaron 
tambien tecnicas de supersimetrfa en el contexto de dinamica biologica de 
microtubules y las propiedades de transporte asociadas a sus paredes de dominio. 
En adicion, se desarrollo una extension supersimetrica compleja del oscilador 
armonico clasico por el cual obtuvimos nuevos modos de oscilacion; resultados que 
pueden extenderse a optica fisica y la fisica de cavidades. Ademas, se desarrollo 
una aplicacion a la nsico-qufmica de moleculas diatomicas usando procedimientos 
de supersimetrfa y de factorizacion. 

La segunda parte contiene un estudio referente a sincronizacion de la dinamica 
caotica de dos neuronas Hodgkin-Huxley, en donde se han aplicado los metodos 
matematicos pertenecientes a la teorfa de control geometrico. Aunque se han 
utilizado diferentes parametros para cada una de las dos neuronas, nuestro estudio 
muestra que se obtienen estados dinamicos de sincronizacion. La sincronizacion 
se logra por la estructura de retroalimentacion de la interconexion (acoplamiento). 
Se muestran los resultados numericos para los estados de dinamica neuronal 
obtenidos. 
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Preface 



Scientific research and technological progress are important characteristics of the 
modern world. They represent fundamental activities that can help mankind to 
understand and transform nature with the puipose of improving standards of life. 
Almost three years have past since I started my doctoral degree activity with the 
hope to contribute myself to the worldwide scientific knowledge. The lines of 
research I chose were on the border between mathematics and biology because I 
was convinced that the interdisciplinary activity is very rewarding and could give 
me better perspectives. 

The doctoral thesis consists of four parts, of which the first contains five 
chapters and is devoted to factorization methods of differential equations and their 
applications in biology and physics, whereas the second part is divided in two 
chapters and deals with the synchronization phenomena as studied in neuronal 
ensembles. The thesis ends up with a final conclusion and the bibliography 
presented in Parts III and IV, respectively. 

The first chapter is a general presentation of the factorization methods for linear 
second order differential equations. Also, the organization for Part I of the thesis 
is presented. 

The second chapter contains an original result for performing factorizations of 
second order differential equations with polynomial nonlinearities that has been 
reported in a paper published in Physical Review E in 2005. At the same time the 
novel procedure allows to obtain particular solutions of travelling kink type in a 
very efficient way. 

The third chapter presents more applications of the method to more complicated 
nonlinear differential equations. The results of this chapter are published in 
Progress of Theoretical Physics in 2005. 

In the fourth chapter, I included the results of a supersymmetric factorization model 
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in the context of microtubules that we published in Physics Letters A in 2003. 

The fifth chapter refers to the original results that have been published in Journal 
of Physics A in December of 2004. A complex extension to the classical harmonic 
oscillator based on a supersymmetric factorization procedure that has been applied 
before in particle physics is introduced in this chapter. The application of the same 
method to the case of Morse potential, a well-known exactly solvable problem 
in quantum mechanics with many applications in the physics and chemistry of 
diatomic molecules is also included here; these results are published in Revista 
Mexicana de Ffsica, 2005. 

With the sixth chapter starts the second part of the thesis. Some remarks on the 
kink type results obtained through factorization methods in the first part for pulse 
propagation along neuron axons, and the connection with the synchronization 
dynamics of a minimal ensemble of two neurons, employing nonlinear control 
theory are presented. 

In the seventh chapter, we focus first on synchronization phenomena from the 
standpoint of their role and importance in natural and technical systems. The 
concept of chaos and the presence of chaotic behavior in nature are also described. 
Next, synchronization methods for the control of chaos and their applications in 
biological systems are shortly reviewed. The problem of the synchronization of two 
Hodgkin-Huxley (HH) neurons is emphasized because of its possible implications 
in the dynamical processes of the brain. A brief discussion of the widely known 
HH mathematical model of the neuron is given. Also, in the Introduction section, 
the organization of Chapters 7 and 8 belonging to Part II of the thesis is presented. 

In the eighth chapter, numerical results for the synchronized dynamics of two 
HH neurons are presented. The mathematical methods employed belong to the 
theory of geometrical nonlinear control and are used with the goal of studying the 
synchronization of two HH neurons that are unidirectionally coupled. These results 
are published in Chaos, Solitons and Fractals in July 2005. 

The order of published papers in this thesis is the following: 

Chapter 2. H.C. Rosu, O. Cornejo-Perez, Supersymmetric pairing of kinks for 
polynomial nonlinearities, Phys. Rev. E 71, 046607 (2005). 

Chapter 3. O. Cornejo-Perez, H.C. Rosu, Nonlinear second order ODE's: factor- 
izations and particular solutions, Prog. Theor. Phys. 114, 533 (2005). 

Chapter 4. H.C. Rosu, J.M. Moran-Mirabal, O. Cornejo, One-parameter nonrel- 
ativistic supersymmetry for microtubules, Phys. Lett. A 310, 353 (2003). 

Chapter 5. H.C. Rosu, O. Cornejo-Perez, R. Lopez-Sandoval, Classical harmonic 
oscillator with Dirac-like parameters and possible applications, J. Phys. A 37, 
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11699 (2004). O. Cornejo-Perez, R. Lopez-Sandoval, H.C. Rosu, Riccati nonher- 
miticity with application to the Morse potential, Rev. Mex. Fis. 51, 316 (2005). 

Chapter 8. O. Cornejo-Perez, R. Femat, Unidirectional synchronization of 
Hodgkin-Huxley neurons, Chaos, Solitons and Fractals 25, 43 (2005). 
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1 



Factorization techniques for linear second order 
differential equations 



1.1 Introduction 



Factorization methods are powerful yet simple algebraic procedures to find 
eigenspectra and eigenfunctions of differential operators that avoid "cumbersome 
transformations, recourse to the ready-made equipment of the mathematical 
warehouse or expansion into power series", to cite from the very first paragraph 
of the 1940's papers of Schrodinger \\\. At the present time, one can find in 
the literature very good informative review papers on the factorization topics 
I2|3l- It is now well known that for second-order linear differential operators, 
the factorizations are equivalent to their Darboux isospectrality (or covariance) and 
also they represent a simple form of intertwining [3]. In this introduction, we will 
touch upon both these issues. 

In the case of Sturm-Liouville operators, E. Schrodinger first developed a 
factorization method he called "that of adjoint first order operators" in 1940- 
1941 Q, during the period he lived in Dublin. In his very first paper on the 
method, Schrodinger deals with four cases: the Planck (harmonic) oscillator, the 
nonrelativistic hydrogen atom, the spherical harmonics in the three-dimensional 
hypersphere, and the Kepler motion in the hypersphere. 
For the quantum harmonic oscillator, he wrote the amplitude equation 

-^--x 2 y + Xy = , (1.1) 
and noticed that it can be written in two different factorized forms 
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Operating on one of these equations with the second of the two first order 
differential operators which occur in it, one gets for the function which results 
from y by applying that operator an equation of the other type, but with A + 2 or 
X — 2, respectively, instead of X. The mutual adjointness of the two first order linear 
operators maintains the quadratic integrability of the solutions and furthermore the 
whole spectrum can be obtained by repeated application of the adjoint operator to 
the X = 1 solutions of the partner operators, e.g., 

^- +x ) ¥ + = y+=e-T (1.2) 



dx 

leads to the odd eigenfunctions in the form 

— -x\ yA+ , (1.3) 

whereas the even eigenfunctions are obtained similarly from the function . 
In his last paper on the method [4|, Schrodinger factorized the hypergeometric 
equation, finding that there are several ways of factorizing it. His factorization 
procedure originated "from a, virtually, well-known treatment of the oscillator", 
i.e., an approach that can be traced back to Dirac's creation and annihilation 
operators for the harmonic oscillator (5) and to older factorization ideas in a paper 
of Pauli |6) and in Weyl's treatment of spherical harmonics with spin Q. It 
should be noted that whereas Dirac's first-order operators were considered only 
as a trick (or 'stratagem'), too insignificant to replace the Sturm-Liouville theory, 
Schrodinger speaks neatly about a method and applies it in a systematic way. 
However, Schrodinger's works were not very much taken into account perhaps 
because of the war years. 

A decade later, in 1951, Infeld and Hull 00 wrote an influential paper in which 
they introduced a different factorization method that became widely known. They 
studied equations of the form 

[M(x,m) + X°]y>;i m (x)=0, 

where M(x,m) is an operator of the form 

d 2 

M(x,m) = -p^ +r(x,m) 

and m = 1,2, ...,n plays the role of a parameter in the potential, whereas the specific 
feature of their method is that the eigenvalue A° is the same for all values of m. 
Infeld and Hull noticed that such equations can be written in two factorized forms 

[- 6+ (m, m + 1)<9_ (m + 1 ,m) - L(m + 1 ) + X%]y%_ m (x) = 

and 

l-d_(m.m- 1)0+ (m- l,m) -L(m) + X^_ m (x) =0 . 
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The eigenfunctions of the neighboring operators M(x,m) and M(jc,m±l) are 
connected by the following relations 

6 + (m-\,m)y':_ m (x) = [A° - L(m + l)]£_ m (*) 

and 

0_(m+l,m)^_ m (x)=j^ +1 (x). 
In addition, the condition 

O_(n + l,n)yS(x) = 

is satisfied leading to 

K=L{n+l). 

The eigenfunctions y®(x) of the operator M(x,0) can be obtained from y^x) 
multiplicatively 

y° n (x)=d + (0,l)6 + (l,2)-d + (n-l,n)f (x) . 

Nothing noteworthy happened for thirty years until Witten |9j wrote a paper 
on dynamical breaking of supersymmetry, in which supersymmetric quantum 
mechanics (SUSYQM) was introduced as a toy model for supersymmetry breaking 
in quantum field theories. 

The SUSY breaking is presented by Witten as a sort of "phase transition" with 
the order parameter being the Witten index, defined as the grading operator 
T = (— 1) f, where Nf is the fermion number operator. For the case of one- 
dimensional SUSYQM, Witten's index operator is the third Pauli matrix (73, which 
is +1 for the bosonic sector and -1 for the fermionic sector of the one dimensional 
quantum problem at hand. It became also quite common to call a particular Riccati 
solution as a (Witten) "superpotential". Papers that now are standard references 
are published during 1982-1984. For example, a breakthrough algebraic result has 
been obtained in 1983 by Gendenshtein [ 10] who introduced the important concept 
of shape invariance (SI) in SUSYQM. The SI property is displayed by some classes 
of potentials with respect to their parameter(s), say a n , and reads 

V n+l (x,a n ) = V n (x,a n+1 )+R(a n ) , 

where R should be a remainder independent of x. This property assures a fully 
algebraic scheme for the spectrum and wave functions. Fixing Eq = 0, the excited 
spectrum is given by the algebraic formula 

n+l 

E„ = £ R{ak) , 

k=2 

and the wave functions are obtained from 

n 

y n (x,ai) = Y[A + (x,a k )y (x,a n+ i) . 
k=\ 
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Another remarkable result of that period is due to Mielnik [11|, who provided 
the first application of the general Riccati solution to the harmonic oscillator, 
obtaining a harmonic potential with an additive tail of the type D 2 [In erf + const.] 
similar to the Abraham-Moses class of isospectral potentials in the area of inverse 
scattering. D. Fernandez gave a second application to the atomic hydrogen 
spectrum, whereas M.M. Nieto clarified further the inverse scattering aspects of 
Mielnik's construction. Mielnik's procedure may be seen as a double Darboux 
transformation in which the general Riccati (superpotential) solution is involved. 
In addition, Andrianov and his collaborators llT2l discovered the relation between 
SUSYQM and Darboux Transformations (DT) or Darboux covariance while 
playing with matrix Hamiltonians in SUSYQM. 



1.2 Darboux covariance 



The Darboux covariance of a Sturm-Liouville equation is clearly stated by 
Matveev and Salle fT3ll . Consider the equation 

and perform the following DT (denoted by i/4l], 

Y^W) = (D-a l )\if=Yx-Oi\if= W ^ ¥ " > , 

¥1 

u — > u[l] =u — 20\ x = u — 2D 2 In \\f\ , 



where 



is the sigma notation of Matveev and Salle for the logarithmic derivative, and W is 
the Wronskian determinant. Then, the Darboux-transformed equation becomes 

-xir xx [l]+u[lMl]=Xy[l] , 

i.e., the spectral parameter X does not change (a result known as Darboux 
isospectrality). When DTs are applied iteratively one gets Crum's result. One 
can also say that the two SL equations are related by a DT. 

Following Matveev and Salle, in order to demonstrate the equivalence of SUSYQM 
with a single DT we consider two Schrodinger equations 

—D 2 \j/ + u\j/ = XY , 
-D 2 + v0 = X<\> , 

related by DT, i.e., v = u[\] and (j) = y[l]> and notice that the function <f>i = 
satisfies the Darboux-transformed equation for X = X\ . 
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If now one uses the second (transformed) equation as initial one and perform the 
DT with the generating function ty\ , one just goes back to the initial u equation. 
That is why one can think of the latter procedure as a sort of inverse DT that can 
be obtained from the direct one as follows: 

M = v-2D 2 ln0i =v[-l] =v-2D 2 ln^f 1 , 

¥ = U ~ ^) (Ai - A) = + ^\ (X 1 - X) . 
Using the sigma notation, 

a _ ¥\x - -tli 
~ ¥1 0i 

the Riccati (SUSYQM) representation of the Darboux pair of Schrodinger 
potentials is obtained 

u = v[-l] = a x + a 2 + Xi , 

v = u[\] = -a x + a 2 + X\ . 

It is now easy to enter the issue of SUSYQM concept of supercharge operators. 
For that, one employs the factorization operators 

B + = -D + o, B =D + o. 

They effect the wave function part of the direct and inverse DT, respectively. 
Moreover, 

B + B = -D 2 + v-Xi , 

B B + = -D 2 + u-Xi . 

Thus, the commutator [B + ,B~]=v — u = —2D 2 In \\f\ gives the Darboux difference 
in the shape of the Darboux-related potentials. Introducing the Hamiltonian 
operators 

H + =B ' B + + Xi , 
H = B + B + X\ , 

one can also interpret the B operators as factorization ones and write the famous 
matrix representation of SUSYQM, as well as the simplest possible superalgebra. 
The factorizing operators in matrix representation are called supercharges in 
SUSYQM, and are nilpotent operators 

fi-=A_a + =( A °_ 0) , (2-) 2 =0, 

and 

Q +=A +°- = {°o A o) ' ( e+ ) 2 = - 
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a = ( jj J J and a + = ( ^ jj J are Pauli matrices. In this realization, the 
matrix form of the Hamiltonian operator reads 



A + A- \ _ / H 
A-A+ " H+ 



defining the partner Hamiltonians as diagonal elements of J^f. They are partners in 
the sense that they are isospectral, apart from the ground state ty gr - of H , which 
is not included in the spectrum of H + . 



1.3 The Mielnik construction 



An interesting possibility to build families of potentials strictly isospectral 
with respect to the initial (bosonic) one arises if one asks for the most 
general superpotential (i.e., the general Riccati solution) such that V + (x) = 
Wg + where V + is the fermionic partner potential. It is easy to see 
that one particular solution to this equation is w p = w(x), where w(x) is the 
common Witten superpotential. One is led to consider the following Riccati 
equation Wg + = wj; + -j^, whose general solution can be written in the form 
Wg(x) = Wp(x) + ^y, where v(x) is an unknown function. Using this ansatz, one 
obtains for the function v(x) the following Bernoulli equation 

■2v(x)w p (x) = l, (1.4) 



dx 



that has the solution 



where J^o(x) = / c x Ug(y) dy, (c = — °° for full line problems and c = for half line 
problems, respectively), and jJ. is an integration constant thereby considered as a 
free parameter. Thus, w g (x) can be written as follows 



d_ 

dx 

w p (x) + a (A) 



w g (x; i u) = w p (x) + — ]n(J r (x)+n) 



d 
dx 



In 



uo(x) 



J^o(x) + M 



(1.6) 
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Finally, one easily gets the V_ (x; ju) family of potentials 



V_(x; M ) = w*(x;m)- b ^' m; 

= V_(x)-2^[ln(^ (x) + At: 

= V_(x)-2(T , x (m) 

v fxl 4u °( x K( x ) i 2u o( x ) n<r> 
1 j S (x)+n + (^o(x)+At) 2 ' ^ ' ^ 

All V_(x;ju) have the same supers ymmetric partner potential V + (x) obtained by 
deleting the ground state. They are asymmetric double-well potentials that may be 
considered as a sort of intermediates between the bosonic potential V_ (x) and the 
fermionic partner V + (x) = V_(x) — 2ob,x(x). From the last rhs of Eq. (11.61) one 
can infer the ground state wave functions for the potentials V_ (x; /J.) as follows 

where f (/x) is a normalization factor that can be shown to be of the form 
f(jli) = \/ju(jU + 1). One can now understand the double Darboux feature of this 
construction by writing the parametric family in terms of their unique "fermionic" 
partner potential 

V_(x;A0=V + (x)-2-^lnf J_ \ (1.9) 



dx 2 \uo(x;ju) 

which shows that the Mielnik transformation is of the inverse Darboux type, 
allowing at the same time a two-step (double Darboux) interpretation, namely, in 
the first step one goes to the fermionic system and in the second step one returns to 
a deformed bosonic system. 

An application of this construction to microtubules is presented in Chapter 4. 



1.4 The connection with intertwining 



Intertwining has been introduced by the French mathematician J. Delsarte 
in 1938 fl4l as an operatorial relationship involving so-called transformation 
(or transmutation) operators but the second World War delayed the detailed 
mathematical studies that came only in the 1950's. By definition, two operators 
Lo and L\ are said to be intertwined by an operator T if 

UT = TU. (1.10) 

If the eigenfunctions <po of Lo are known, then from the intertwining relation one 
can show that the (unnormalized) eigenfunctions of L\ are given by (pi = T(pQ. The 
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main problem in the intertwining transformations is to construct the transformation 
operator T . One-dimensional quantum mechanics is one of the simplest examples 
of intertwining relations since Witten's transformation operator T qm = 7i is just a 
first spatial derivative plus a differentiable coordinate function (the superpotential) 
that should be a logarithmic derivative of the true bosonic zero mode (if it exists), 
but of course higher-order transformation operators can be constructed without 
much difficulty. 

Thus, within the realm of the one-dimensional quantum mechanics, writing T\ = 
D — — , where u is a true bosonic zero mode, one can infer that the adjoint operator 

7y = —D — ^ intertwines in the opposite direction, taking solutions of L\ to those 
of L 

<ft) = 7\ + <pi . (1.11) 

In particular, for standard one-dimensional quantum mechanics, Lo = H and 
L\ = H + and although the true zero mode of H is annihilated by T\, the cor- 
responding (unnormalized) eigenfunction of H + can nevertheless be obtained by 
applying T\ to the other independent zero energy solution of H . It is only in 
the last decade or so, that the intertwining approach becomes well-known to the 
SUSYQM factorization community and some authors start to play with higher- 
order generalizations. But, as always, the most important (at least for standard 
quantum mechanics) are the simplest cases, namely the Darboux first-order inter- 
twining operators. 

The first part of this thesis deals with factorization methods, among which an 
original factorization of nonlinear second order ordinary differential equations 
(ODE) and supersymmetric techniques, as applied to some biological and physical 
systems. Chapters 2 and 3 contain explicitly the new factorization procedure 
developed by us to obtain kink type solutions for nonolinear second order ODE 
that describe several important processes, for instance, the tubulin polymerization 
in microgravity conditions and the pulse propagation along nerve axons. In Chapter 
4, supersymmetric approaches are applied in the framework of biological dynamics 
of microtubules (MTs); the latter results are related to transport properties 
associated to the MT domain walls. In Chapter 5, applications of supersymmetric 
factorization procedures in some physical systems are presented. A complex 
extension for the classical harmonic oscillator by means of a direct relationship 
between the Dirac and Schroedinger equations is obtained. In addition, the 
same procedure is applied to a molecular physics problem in connection with the 
dissociation of diatomic molecules. 
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2 



A new factorization technique for differential equations 
with polynomial nonlinearity 



Abstract. In this chapter, it is shown how one can obtain kink solutions of ordinary 
differential equations with polynomial nonlinearities by an efficient factorization procedure 
directly related to the factorization of their nonlinear polynomial part. This is different 
of previous factorization procedures of differential equations of this type that have been 
performed by only a few authors, most notably by Berkovich 1171 . Of main interest 
here because of their numerous applications are the reaction-diffusion equations in the 
travelling frame and the damped-anharmonic-oscillator equations. In addition, interesting 
pairing of the kink solutions, a result obtained by reversing the factorization brackets in 
the supersymmetric quantum mechanical style, are reported. In this way, one gets ordinary 
differential equations with a different polynomial nonlinearity possessing kink solutions of 
different width but propagating at the same velocity as the kinks of the original equation. 
This pairing of kinks could have many applications. The mathematical procedure is 
illustrated with several important cases, among which the generalized Fisher equation, the 
FitzHugh-Nagumo equation, and the polymerization fronts of microtubules (MTs). In the 
latter case, a new polymerization front is predicted that can show up in solutions containing 
MTs borne on satellites. Because of the microgravity conditions the polymerization rates 
could deviate from the normal ones and this could lead to a change of the width of the 
polymerization front. 



2.1 Introduction 



Factorization of second-order linear differential equations, such as the Schrodinger 
equation, is a well established method to get solutions in an algebraic manner 
EIU El- We are interested in factorizations of ordinary differential equations 
(ODE) of the type 

u" + yu'+F(u) = , (2.1) 
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where F(u) is a given polynomial in u. If the independent variable is the 
time then 7 is a damping constant and we are in the case of nonlinear damped 
oscillator equations. Many examples of this type are collected in the Appendix 
of a paper of Tuszyhski et al. fToll . However, the coefficient 7 can also play 
the role of the constant velocity of a travelling front if the independent variable 
is a travelling coordinate used to reduce a reaction-diffusion (RD) equation to 
the ordinary differential form as briefly sketched in the following. These RD 
travelling fronts or kinks are important objects in low dimensional nonlinear 
phenomenology describing topologically-switched configurations in many areas 
of biology, ecology, chemistry and physics. 
Consider a scalar RD equation for u(x,t) 

(2.2, 

where Q is the diffusion constant and s is the strength of the reaction process. 
Eq (12.21) can be rewritten as 



du_^_ 
dt dx 



+ F(u) , (2.3) 



where the coefficients have been eliminated by the rescalings t = st and x = 
{s/@y' 2 x, and dropping the tilde. Usually, the scalar RD equation possesses 
travelling wave solutions w(£ ) with E, = x — vt, propagating at speed v. For this 
type of solutions the RD equation turns into the ODE 

u" + vu'+F(u) = , (2.4) 

where ' = D = j?. Eq. (12.4b has the same form as nonlinear damped oscillator 
equations with the velocity playing the role of the friction constant. 
For applications in physical optics and acoustics it is convenient to write the 
travelling coordinate in the form t, = kx — cot = k(x — vt) with k\ = 00. This is 
a simple scaling by k of the previous coordinate turning Eq. (12.41) into the form 

u" + \' + ^F(u) =0 (2.5) 
k k 

that can be changed back to the form of Eq. ( 12. II ) by redefining f = j and 
F(u) = ^F(u). 

In general, performing the factorization of Eq. (12.11) means the following 



D-h(u) D — fi (u) 



. (2.6) 



This leads to the equation 



u" p-uu! — f\u — fovl + fxfiu = . (2.7) 

du 
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The following groupings of terms are possible related to different factorizations: 
a) Berkovich grouping: In 1992, Berkovich jTTl proposed to group the terms as 
follows 

u" - (fi + h) u + L/i/2 - u = 0, (2.8) 

and furthermore discussed a theorem according to which any factorization of an 
ODE of the form given in Eq. (12.61) allows to find a class of solutions that can be 
obtained from solving the first-order differential equation 

u'-f x {u)u = Q. (2.9) 

Substituting the first-order ODE (I2.9t in the Berkovich grouping one gets 

(fib + fib) u' + (fufu, - ^f/i*") u = 0, (2. 10) 



u" 



where we redefined f\ (u) = f\b{u) and fi(u) = fib(u) to distinguish this case from 
our proposal following next. For the specific form of the ODEs we consider here, 
Berkovich's conditions read 

fib(-7-fib-^u)=^, (2.11) 
V du J u 



fYb + fv, = -y. (2.12) 
b) Grouping of this work: We propose here the different grouping of terms 

' ^P-U + 01+^2 )«' + 01^2" = (2.13) 



\ du 

that can be considered the result of changing the Berkovich factorization by setting 
fib(u) = 0i (u) and fzb{u) — > <h.{ u ) under the conditions 

F(u) 

0102 = ^, (2.14) 



0i +02 + ^u = -y. (2.15) 
du 

The following simple relationship exists between the factoring functions: 

02 W =hb\u) u 

du 

and further (third, and so forth) factorizations can be obtained through linear 
combinations of the functions fib, f2b and 02- 

Based on our experience, we think that the grouping we propose is more 
convenient than that of Berkovich and also of other people employing more 



16 



difficult procedures. The main advantage resides in the fact that whereas in 
Berkovich's scheme Eq. ( 12.1 It is still a differential equation to be solved, in 
our scheme we make a choice of the factorization functions by merely factoring 
polynomial expressions according to Eq. d2.14t and then imposing Eq. ( I2.15t leads 
easily to an ^-depending 7 coefficient for which the factorization works. This fact 
makes our approach extremely efficient in finding particular solutions of the kink 
type as one can see in the following. 

In the next section, it is shown on the explicit case of the generalized Fisher 
equation all the mathematical constructions related to the factorization brackets and 
their supersymmetric quantum mechanical like reverse factorization. In less detail, 
but within the same approach, damped nonlinear oscillators of Dixon-Tuszyriski- 
Otwinowski type and the FitzHugh-Nagumo equation, are studied in Sections 2.3 
and 2.4, respectively. 



2.2 Generalized Fisher equation 



Let us consider the generalized Fisher equation given by 

u" + yu +u(l-u n ) =0, (2.16) 

The case n = 1 refers to the common Fisher equation and it will be shortly 
discussed as a subcase. Eq. i2. 14b allows to factorize the polynomial function 



Now, by choosing 



{\-u n ) = {\-u n ' z ){\ + u n l l ), 



<j> 1 =a 1 (l-u n / 2 ), 4> 2 = 1(1 + M »/ 2 ) , ai^O, 

the explicit forms of a\ and 7 can be obtained from Eq. d2. 15l > 

^l u + ^ + ^ 2 = -^a l u n l 2 +a x (l-u n l 2 ) + (l/a l ){\+u n l 2 ) = -y. 
au 2 

Introducing the notation h n = (§ + l) 1 ^ 2 one gets: a\ = ±h~ l , 7 = =F {K + K 1 ) ■ 
Then Eq. (12.161) becomes 



u" ± {hn + h~ l ) u' + «( 1 - u") = 
and the corresponding factorization is 



D ± h n {u n l 2 + 1)1 \d T h-\u n l 2 - 1) 



w = . 



(2.17) 



(2.18) 



It follows that Eq. d2.17t is compatible with the first-order differential equation 

u'^K 1 (u"/ 2 -l)u = . (2.19) 
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Integration of Eq. ( 12.191 ) gives for 7 > 

'l±e X p[(/i, I -/ I „- 1 )(^-^ ) 
Rewritten in the hyperbolic form, we get 
1 1 



-2/n 



(2.20) 



■tanh 

2 2 

1 1 . 
- — - coth 

2 2 



L2 



I n \ 2/n 



1 {K-h-') (£-§,)" 



(2.21) 



The tanh(-) form is precisely the solution obtained long ago by Wang fTE\ and 
Hereman and Takaoka [ 19 1 by more complicated means. 
Moreover, a different solution is possible for 7 < 



or 



1 ± exp 



1 1 . 

2 + 2 tanh 

1 1 . 

2 + 2 C ° th 



(h n -h- l )(^-^) 



-2/n 



(2.22) 



1 ,\V" 

.-(hn-h- 1 )^-^)^ 



(2.23) 



respectively. 



2.2.1 Reversion of factorization brackets without the change of the scaling 
factors 

Choosing now 0i = a\(l + u n l 2 ) and 02 = jjj- (iW 2 — l) leads to the same equation 
(12.171) but now with the factorization 



D^h n (u n/2 -l) D±h-\u n/2 + \) u = 0, 
and therefore the compatibility is with the different first-order equation 

u'±K l (u n/2 + l)u = . 



(2.24) 



(2.25) 



However, the direct integration gives the solution (for 7 > 0) 

/ \ 2 l n 



\ \±tX V [{h n -hn l )^-^)]j 

(-If/" (l± e x V [{h n -h- l )(^-^ 



-2/n 



(2.26) 



which are similar to the known solution Eq. (I2.20I) . For 7 < 0, solutions of the type 
given by Eq. ( I2.22t are obtained. 
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2.2.2 Direct reversion of factorization brackets 

Let us perform now a direct inversion of the factorization brackets in d2.18t similar 
to what is done in supersymmetric quantum mechanics in order to enlarge the class 
of exactly solvable quantum hamiltonians 



DtK\u^ 2 -\)] \D±h n (u"/ 2 + \ 



. 



(2.27) 



Doing the product of differential operators the following RD equation is obtained 



n/2 



\-hiu n ' 2 



u" ± (h n + h n ) i/ + w 1 + u 
Eq. d2.28t is compatible with the equation 

u'±h n (u n / 2 + i)u = 0, 
and integration of the latter gives the kink solution of Eq. ( 12.281 1 



0. 



(2.28) 



(2.29) 



1 



l±expp3-fc n )(£-& 



1 ±exp 



(fi-hn)(S-Zo) 



(2.30) 

for 7 > 0. On the other hand, for y < the exponent is the same but of opposite 
sign. Hyperbolic forms of the latter solutions are easy to write down and are similar 
up to widths to Eqs. (12.211) and (I2.23I) . respectively. 

Thus, a different RD equation given by (12.281) with modified polynomial terms and 
its solution have been found by reverting the factorization terms of Eq. 



Although the reaction polynomial is different the velocity parameter remains the 
same. The main result, which is a general one, that we find here is the following: 
At the velocity corresponding to the travelling kink of a given RD equation there 
is another propagating kink corresponding to a different RD equation that is 
related to the original one by reverse factorization. We can call this kink as the 
supersymmetric (susy) kink because of the mathematical construction. 
Finally, one can ask if the process of reverse factorization can be continued with 
Eq. (I2.28I) . It can be shown that this is not the case because Eq. (12.281) has already a 
discretized (polynomial-order-dependent) y and this fact prevents further solutions 
of this type. Suppose we consider the following factorization functions 



\-hiu n ' 2 



02 = 3\ ( 1 + u 



(2.31) 



Then, one gets d\ = and solve a^ 1 +d\ = h~ l +h n . The solutions are: n = 0, 
which implies linearity, and n = —4, which leads to a Milne-Pinney equation. On 
the other hand, Eq. ( 12.281 ) with an arbitrary y can be treated by the inverse factor- 
ization procedure to get the susy partner RD equation and its susy kink. 
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2.2.3 Subcase n = 1 

This subcase is the original Fisher equation describing the propagation of mutant 
genes 

du d 2 - 



- = — + u(l-u). 
In the travelling frame, the Fisher equation has the form 

u" + yu' + u(l - u) =0 . 



(2.32) 



(2.33) 



When the y parameter takes the value Yi — f \/6 (i.e., h\ = ^) one can factor 
Fisher's equation and employing our method leads easily to the known kink 
solution 

1 



Up = - I 1 — tanh 



(2.34) 



that was first obtained by Ablowitz and Zeppetella fTTl with a series solution 
method. On the other hand, the susy kink for this case reads 



susy 



1 

4 



1 - tanh 



V6 



(2.35) 



i.e., it has a width one and a half times greater than the common Fisher kink and is 
a solution of the partner equation 



// . 5\/6 , 
u H — 



5 A/2-l 







(2.36) 



A plot of the kinks up and «F,susy is displayed in Fig. 2. 1 




Fig. 2.1: The front of mutant genes (Fisher's wave of advance) in a population and the 
partner susy kink propagating with the same velocity. The axis are in arbitrary units. 
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2.2.4 Subcase n = 6 

This subcase is of interest in the light of experiments on polymerization patterns 
of MTs in centrifuges. It has been discovered that the polymerization of the 
tubulin dimers proceeds in a kink-switching fashion propagating with a constant 
velocity within the sample. Portet, Tuszynski and Dixon |20j used RD equations 
to discuss the modification of self-organization patterns of MTs as well as the 
tubulin polymerization under the influence of reduced gravitational fields. They 
used the value n = 6 for the mean critical number of tubulin dimers at which the 
polymerization process starts and showed that the same nucleation number enters 
the polynomial term of the RD process for the number concentration c of tubulin 
dimers 



c" + -c' + c l 



1 



5 )=0 



(2.37) 



The polymerization kink in their work reads 



cptd 



1 -tanh 



1/3 



(2.38) 



On the other hand, the susy polymerization kink (see Fig. (2.2)) of the form 



■'Susy 



2- l-tanh[3($-$ )] 



1/3 



can be taken into account according to the hyperbolic form of Eq. 
propagates with the same speed and corresponds to the equation 




c"±-c' + c 



1- 15c 3 -16c 6 ) = 0. 



(2.40) 



In principle, this equation could be obtained as a consequence of modifying the 
kinetics steps in the microtubule polymerization process. 



Tubulin \ \ 










\/ CpTD 


V 

- wave front 


\ 

\ c 

1 \ / susy 







-5 5 10 15 



Fig. 2.2: The polymerization kink of Portet, Tuszynski and Dixon |i20J and the susy 
kink propagating with the same velocity. 
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2.3 Equations of the Dixon-Tuszynski-Otwinowski type 



In the context of damped anharmonic oscillators, Dixon et al. l2"2l studied 
equations of the type (in this section, we use ' = D T = 4-) 



it . I , a n—1 a , i 

U +U +AU — U =U -\-U +U 



(VA-u L i- l )(VA + u?- l )=0 (2.41) 



and gave solutions for the cases A = | and A = with n = 4 and n = 6, 
respectively. For this case, time is the independent variable. The factorization 
method works nicely if one uses g n = \fnj2 and dealing with the more general 
equation 

u"±VA(g n + g- l )u' + u(A-u n - 2 ) =0, (2.42) 
for which we can employ either the factorization functions 



or 



1-1 



Then, Eq. (12.421) can be factored in the forms 

D t ±g n (u'i- l + VA)] \D^g- l (u L 2- 1 -VA) 



and 



D t Tgn{u'i- l -VA)] \D z ±g-\ui- l + ^A) 



u = 



u = 



Thus, Eq. <!2.42b is compatible with the equations 

"Tg,; 1 {u^ 1 - vXj u=o, 

u'±g~ l (u'i- 1 + VX)u = 



(2.43) 
(2.44) 

(2.45) 
(2.46) 



that follows from Eq. (12.43b and Eq. (12.44b . Integration of Eqs. ( 12.45b . (12.46b gives 
the solution of Eq. (12.42b 



Va 



1 ±exp 



and 



VA(g n -g„ l )(z-ZQ 



Va 




7>0 



(2.47) 



1 ±exp 



VA(g n -g n l )(T-T Q 




7<0 . 



(2.48) 
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The solutions obtained by Dixon et al. are particular cases of the latter formulas. 
Reversing now the factorization brackets in (12.431) 



fiT?„' (u'i- l -VA)\\D t ±g n (u"2- l + VA) 
leads to the following equation 

u" ± VA(g n +gn l )u + U (VX + M? -1 ) ^VA-^-U^ 1 

which is compatible with the equation 

whose integration gives the solution of Eq. ( 12.501 ) 



u = 



(2.49) 



: 0, (2.50) 



(2.51) 



and 



\l±exp[VAg n (x-Xo)]^ 

Va 

l±exp[-y/Ag n (x-Zo)] 



i 

n-2 



7>0 



7<0. 



(2.52) 



(2.53) 



2.4 FitzHugh-Nagumo equation 



Let us consider the FitzHugh-Nagumo equation, which is a common approxima- 
tion to describe nerve fiber propagation, 



du d u 
dt dx 2 



+ h(1 — u)(a — u) = , 



(2.54) 



where a is a real constant. Moreover, if a = — 1, one gets the real Newell- 
Whitehead equation describing the dynamical behavior near the bifurcation point 
for the Rayleigh-Benard convection of binary fluid mixtures. The travelling frame 
form of (12.541) has been discussed in detail by Hereman and Takaoka ff9l 

u" + yu' + u(u — I) (a — u) = 0. (2.55) 

The FitzHugh-Nagumo polynomial function allows the following factorizations: 

0! =±{V2)- 1 (u-\), <fo = ±Vl{t 



[a — u) 



when the y parameter is equal to y a \ = ± ^ that we also write as y a ^\ 
±y/a{ga\-ga\), where gal = 



2a. 
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In addition, we can employ the factorization functions 

01 = ±{V2y l (a-u), fc = ±V2(«-l) 

when 7 a 2 = db ~*^ , or written again in the more symmetric form y Q) 2 = ±v^(#a2 — 
\fajl. Thus, Eq. ( I2.55t can be factored in the two cases 

u" ± y a ,\u' + w(w — — u) = 0, (2.56) 



g^ 1 ), where g a2 



and 



w" ± y a ,2U + w(w — l)(a — u) = . 



(2.57) 



In passing, we notice that for the Newell -Whitehead case a = — 1 the two equations 
coincide and are the same as the generalized Fisher equation for n = 2. 
In factorization bracket forms, Eqs. d2.56t and ( 12.571 ) are written as follows 



£>=F v2(a-u) 



D±(V2)- l (\ 



u = 



and 



DT\/2(m-1) D=f(V2y\a-u) u = 0, 
and are compatible with the first order differential equations 

u ±{Viy\l-u)u = 0, for Yal , 



i't(V2) 1 (a-u)u = 0, for 7 a , 2 • 



(2.58) 
(2.59) 

(2.60) 
(2.61) 



Integration of Eqs. (I2.60t and (I2.6U gives the solution of Eq. (12.55b for the two 
different values of the wave front velocity y a \ and Ya,2- 
For Eq. (12.561) we get 



1 



1 ±exp[(V2)-i($ - &)] ' 1 ±exp[-(V2)-i (£ - £„)] ' 

for 7„ i positive and negative, respectively. 



(2.62) 



As for Eq. (12.57b . the solutions are 

a 



(2.63) 



1 ±exp[-(v/2)-y§ - So)] ' 1 iexpKv^)- 1 ^ " M ' 

for 7a,2 positive and negative, respectively. 

Considering now the factorizations (12.581) and (12.59b . the change of order of the 
factorization brackets gives 



and 



D±(V5) _1 (1 -u 
DT(V2y\a-u) 



D^v2(a — u 
D=fV2(«-1) 



u = 



u = . 



(2.64) 
(2.65) 



24 



Doing the product of differential operators (and considering the factorization term 
u' — 02" = 0) gives the following RD equations 



and 



u" ± y a \u + u(4u — l)(a — u) = 0, 



u" ± Yaiu' + u(u — l)(a — u — 3u 2 ) = 0, 



Eqs. (I2.66t and (12.67b are compatible with the equations 

u' =p \fl{a — u)u = 

and 

u' ^ \[7.{u — \)u = , 



(2.66) 
(2.67) 

(2.68) 
(2.69) 



respectively. Integrations of Eqs. ( 12.681) and ( 12.69b give the solutions of Eqs. (12.66b 



and (I2.67I) . respectively. The explicit forms are the following: 
(i) for d2~66Ti 

a a 

u > - 



l±exp[-v/2a(€-€o)] 



(ii) for (12371 

M> = 



l±exp[V2(^-<§o)] 



l±exp[V2^-£ )] 



l±«p[-^-&)] 



(2.70) 



(2.71) 



2.5 Conclusion of the chapter 



In this chapter, we have been concerned with stating an efficient factorization 
scheme of ODE with polynomial nonlinearities that leads to an easy finding of 
analytical solutions of the kink type that previously have been obtained by far 
more cumbersome procedures. The main result is an interesting pairing between 
equations with different polynomial nonlinearities, which is obtained by applying 
the susy quantum mechanical reverse factorization. The kinks of the two nonlinear 
equations are of different widths but they propagate at the same velocity, or if 
we deal with damped polynomial nonlinear oscillators the two kink solutions 
correspond to the same friction coefficient. Several important cases, such as 
the generalized Fisher and the FitzHugh-Nagumo equations, have been shown to 
be simple mathematical exercises for this factorization technique. The physical 
prediction is that for commonly occurring propagating fronts, there are two kink 
fronts of different widths at a given propagating velocity. Moreover, the reverse 
factorization procedure can be also applied to the Berkovich scheme with similar 
results. It will be interesting to apply the approach of this work to the discrete 
case in which various exact results have been obtained in recent years l23ll . More 
general cases in which the coefficient y is an arbitrary function are also of much 
interest because of possible applications. The same factorization scheme as it 
works for more complicated ordinary differential equations is described in the next 
chapter. 
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3 



Application to more general nonlinear differential 
equations 



Abstract. In the previous chapter we considered the coefficient in front of the first 
derivative as a constant quantity. However, the employed factorization technique can be 
used almost unchanged for the more general case when the condition of constancy of this 
coefficient is relaxed. In this chapter, we obtain more kink type solutions through the 
same factorization procedure for a number of more general nonlinear ordinary second order 
differential equations with important applications in biology and physics. 



3.1 1 Introduction 

Considering the following type of differential equation 

u" + g(u)u +F(u)=0 (3.1) 

where again as in the previous chapter ' means the derivative D = £ and S, = x — vt ; 
one can factorize Eq. d3.lt in the following form 

[D-fc(u)][D-(j) 1 (u.)}u = 0. (3.2) 

Performing now the product of differential operators leads to the equation 

n d(j){ III 

U — UU — 01 M — 02M +0102^ = 0, (3.3) 

du 

for which one way of grouping the terms is as follows 

u" - ^01 + 02 + ^""J u ' + 01 <h u = • ( 3 - 4 ) 
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Eqs. (I3.lt and (13.41 are lead to the conditions 

g(u) = -(<j> l + <h + ^u) (3.5) 
\ au J 

and 

F '(«) = ^>i02« • (3-6) 

If is a polynomial function, then g(u) will have the same order as the bigger 
of the factorizing functions 0i(w) and tyi{u), and will also be a function of the 
constant parameters provided by the function F(u). 

In the context of classical mechanics, Eq. d3.lt could be seen as an anharmonic 
oscillator with nonlinear damping. The case g(u) = V where v is a constant 
value has been presented in the previous chapter. There, by means of a 
simple factorization method exact particular solutions of the kink type for 
reaction-diffusion equations and damped-anharmonic oscillators with polynomial 
nonlinearities have been obtained. In addition, SUSYQM-like reversing of 
factorization brackets has been performed providing new kink solutions for 
equations with different polynomial nonlinearities. 

Based on the given grouping in Eq. (13.41) for Eq. (I3.lt . a simple mathematical 
procedure is proposed by which one gets particular solutions through factorization 
methods that allows finding solutions satisfying a compatible (nonlinear) first order 
differential equation. 

The puipose of this chapter is to further apply this mathematical scheme to a wealth 
of important cases for which explicit particular solutions are not easy to find in the 
literature or are obtained by more involved techniques. The examples we present 
herein are the modified Emden equation, the Generalized Lienard equation, the 
convective Fisher equation, the generalized Burgers-Huxley equation, all of whom 
have significant applications in nonlinear physics. Explicit particular solutions are 
presented. 



3.2 1 Modified Emden equation 

Let us consider the following modified Emden equation 

u" + auu' + jSw 3 = . (3.7) 
The polynomial F(u) = fiu 3 allows the following factorizing functions 

01 (w) =a\\f]iu , and §2(11) = —\ffiu , 

where a\ is an arbitrary constant. Eq. (13.51) is used to obtain the function g(u), 

g(u)=-(2a iy /pu+±-y/pi?J =-y/p(??±±±^u , 
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then identifying a = — a//3 (j-^—^J (or ci\ + = )» where we use #1 as 

a fitting parameter providing that ai < for a > 0. We note that = g(p,a\;u). 
Eq. (13.71) is now rewritten in the following form 

M "_ v /]3^±iW + j 3 M 3 = 0) (38) 

the equation can be factorized as follows 

\D — ^-\f]5u \ (l) — ai yfpuj u = (3.9) 

and therefore the compatible first order differential equation is 

u - ai y/pu 2 = . (3.10) 
Integration of Eq. ( 13.101 ) gives the particular solution of Eq. d3.8t 

u= = , (3.11) 

where ^ is an integration constant. If we consider the quadratic equation for a\, 
then Eq. (13. Ill is expressed as a function of a, 

4 



(3.12) 




Fig. 3.1: Real part for the factorization curve of the parameter a i =ai + (a,/J) that 
allows the factorization of Eq. Ol . a^O. ae [-10, 10] and /3 G [-10, 10]. 
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Fig. 3.2: Imaginary part for the factorization curve of the parameter a\ + = a\ + (a,j5) that 
allows the factorization of Eq. Ol . fl^O.ae [-10, 10] and j3 g [-10, 10]. 



Let us consider now another pair of factorizing functions 

(l>l(u) =a l yf]3u 2 , (j> 2 (u) = — y/]3 , 



a i 



then, using Eq. (13.5b . the function g(u) = — \f$ M^- + 3aiu 2 j is easily obtained. 
Therefore, the original modified Emden equation (13.71) becomes 

u-\f$ (J^ + 3aiu 2 ^ju' + fiu 3 = . (3.13) 

This equation allows the factorization 

{D-^-yfp) (D-a iy fpu 2 ^u = , (3.14) 

where from we obtain the compatible first order differential equation 

u' -axy/pu 3 = (3.15) 

with the solution 

(3.16) 



The above example shows that different factorizations of F{u) would yield differ- 
ent forms for the function g(u). This is an important consequence of applying this 
mathematical technique to the case g(u) ^ const. 
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3.3 Generalized Lienard equation 



Let us consider now the following generalized Lienard equation with a cubic 
polynomial function F(u) 



u" + g{u)u' + Au + Bu 2 + Cu i = 



(3.17) 



The polynomial function F(u) can be factorized in several ways, we consider first 
the factorization F(u) = u(a + b + Cu) (d — e + u) where 

a = B/2, b = \ / B 2 -4AC/2, d = B/2C, e = \Jb 2 - 4AC/2C, 
and the condition B 2 — A AC > holds. If we consider the factorizing functions as 

01 (u) = a\ (a + b + Cu) and 02 ( M ) = — (d — e + u), 



a\(a+b)+(d-e) ( 



where again a\ is an arbitrary constant that can be used as a fitting parameter, the 
function g(u) 



will be obtained. Then, Eq. (I3.17t is 

u' +Au + Bu 2 + Cu 3 =0 , (3.18) 



rewritten as 

,2 



U + 



a\(a + b) + (d-e) (2a\C+\ 



Cl\ 



+ 



and the corresponding factorization will be 
1 



D- —(d-e + u) 



[D — a\ (a + b + Cu)] u = 



where from the compatible first order differential equation is obtained 

u — a\(a + b + Cu)u = , 

and whose solution is 

= {a + b)exy[ai(a + b){E, -g )] 
1-Cexp[ai (a + 6) 



(3.19) 



(3.20) 



(3.21) 



where (a + b) = B+ ^ B ^ 4AC . 

Let us consider now the following reduction of terms in Eq. (I3.17I) . B = and C = 1 
in order to calculate a particular solution for the so-called autonomous Duffing- van 
der Pol oscillator equation J25j, 



u" + (G + Eu 2 )u' +Au + u 3 = , 



(3.22) 



where G and E are arbitrary constant parameters. The polynomial function allows 
the following factorizing functions 



fi(u) = a\(A + u 2 ) and 02 (m) 



1 
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then g(u) = — y- L ^ l 1" ~iaiu 2 j . Eq. d3.22t is now rewritten 

„ / a 2 A + 1 „ ,\ , , 
u - — h 3a\u l u +Au + = . 

The corresponding factorization of Eq. d3.23t is given as follows 



1 

D 

a\ 



[D-ai(A + u 2 )] u = 



and the obtained compatible first order equation 

u' — a\ (A + u 2 )u = . 



Integration of Eq. d3.25t gives the particular solution of Eq. (13.231 



r ( exp[2aiA(g - gp)] 
U V \,l-exp[2fliA(§-&) 



i/2 



(3.23) 



(3.24) 



(3.25) 



(3.26) 



Comparing Eqs. <!3.22b and (13.23b . «i = — j and G 
( 13.261 ) is now written as a function of A and E, 



AE A +9 
3E 



are obtained. Solution 



u = ±VaC ex P[-i A£ (£-&>)] 



1/2 



l-exp[-|A£(§-§,)] 



(3.27) 



This is a more general result for the particular solution than that obtained by Chan- 
drasekar et al. in 1231 by other means, in fact, it is recuperated when E = j8 and 

A — A- 
p 2 ' 




Fig. 3.3: Real part for the factorization curve of the parameter E + = E + (G,A) that allows 

Note that a \ = 

Ae [-10,10]. 



the factorization of Eq. 13.231 . Note that a\ = -|; E ^ 0. G G [-10, 10] and 
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G 



Fig. 3.4: Imaginary part for the factorization curve of the parameter E + = E + (G,A) that 
allows the factorization of Eq. I ET231 . E^Q.Ge [-10, 10] and A E [-10, 10]. 



3.4 Convective Fisher equation 



Let us consider the convective Fisher equation given in the following form [26 1, 



Id 2 



du 



dt 2dx^ +u{1 - u) -^ U dx 



(3.28) 



where /J. is a positive parameter that serves to tune the relative strength of 
convection. If the variable transformation B, = x — vt is performed, then we obtain 
the following ordinary differential equation 



u' + 2{y - liu)u +2u(\ - u) = . 
The polynomial function allows the factorizing functions 

01 (m) = v / 2a , i (1 — u) and §i{u) = — , 



(3.29) 



and Eq. (I3.5t gives the function g{u 
rewritten as follows 

,2 



-y/l (-^— — 2a\u). Eq. (I3.29t is 



" + 2 f - + V^iM^) w' +2u(1 - «) = . (3.30) 



If we set the fitting parameter a\ 
factorized in the following form 



then we obtain v = ^i 2 .. Eq. d3~30l is 



V2 



D 



V2 



D-V2ai(l 



u = 



(3.31) 
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that provides the compatible first order equation 

u — Vla\u(\ — u) = u + }iu(\ — u) = (3.32) 
whose integration gives 

(3.33) 



l±expLu(£-So)] 



14 ■ 
12 ■ 
1C 



3 
> 



2 4 6 8 10 12 14 



Fig. 3.5: Factorization curve of the parameter v = v(jli) that allows the factorization of 



Eq. ll3~30l . a x = 



3.5 Generalized Burgers-Huxley equation 



In this section we obtain particular solutions for the generalized Burgers-Huxley 
equation discussed by Wang et al. in \ 21\ 



du s du d 2 u s s 

at ox ox 2 - 



(3.34) 



If the coordinates transformation £, = x — vt is performed then Eq. (I3.34t is 
rewritten in the following form 



u" + (v + au b )u + pu(l - u d )(u d - y) = , 
and the polynomial function allows the choice for the factorizing terms 



(3.35) 



0i (m) = v)3ai(l — u s ) and tyi{u) 



(u°-y). 



Eq. (I3.5t provides g(u) = \[$ (j-^- 1 + "'^^ - u s ^\ , and we can do the following 



identification of constant parameters 



a 



-a\{\ + 8) + \ 
a\ 
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Writing Eq. (I3.35t in factorized form 



-^-(u s -r) 



D— \/j8fli(l — u h 



u = , 



the solution 



1/5 



Vl±exp[-flivW§-&)], 
of the compatible first order equation 

u — y/~[}aiu(l — u S ) = , 



(3.36) 



(3.37) 



(3.38) 



is also a particular kink solution of Eq. (I3.35t . Solving the quadratic equation 
( 13.361 ) for c?! = ci\ (a, j3, 8) we obtain 



a\,_ 



-a± V« 2 + 4j3(l + 5) 



2Vi8(l + 5) 

then Eq. ( 13.371 ) becomes a function m = u(oc,P,8;t), and v = v(a,/3,y, 5). If we 
set 5 = 1 in Eq. (13.34b . then we obtain the following particular Burgers-Huxley 
solution 

u = — = , (3.39) 

l±exp[- fll VF(^-§o)] 



andai + _ = v H , v = v(a,j3,y). 




Fig. 3.6: Real part for the factorization curve of the parameter «i = #i + (a, j3, 5 = 1) that 
allows factorization of Eq. J3~35t with g = 1 . a x ^ 0. a 6 [-20,20] and J3 G [-20,20]. 
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Fig. 3.7: Imaginary part for the factorization curve of the parameter 
a \ + = fli + (o!,j3,5 = 1) that allows factorization of Eq. (I3.35> with S = 1. a\ ^ 0. 

a e [-20,20] and j3 e [-20,20]. 



If we chose now the factorizing terms as 

(j>i(u) = \f]3ei(u s — y) and <fc(w) 



(l-u 5 ) 



we obtain = y / /3 " ^'^ 1 + 1 ^"^ m 5 ^, and the following identification 

Eq. <I3.35I > is then 



of parameters v 



and a 



l-ef(l+5) 



factorized in the different form 



D 



D-y^ ei (u S -Y) 



u = 



The corresponding compatible first order equation is now 

u - \ffie\u(u s - y) = , 



(3.40) 



(3.41) 



and its integration gives a different particular solution for Eq. <I3.35I > from that 
obtained for the first choice of factorizing terms d3.36t . however, we point out 
that the parameter a has changed for the second choice of factorizing terms. The 
solution of Eq. ( 13.411 ) is given as follows 



( 7 \ V5 

U= ^liexpfieiViWS -£<>)]/ 

Solving the quadratic equation for e\ = <?i(a,j8,<5) we obtain 



(3.42) 



a± yja 2 + 4/3(1 + 8) 
VF(T+S) 
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then Eq. (13.42b becomes u = u(a,p , y, <5;t), and v = v(a,j3,y, 8). If we set 8 = 1 
in Eq. 03.341 1. then the following Burgers-Huxley solution is obtained 




(3.43) 



Fig. 3.8: Real part for the factorization curve of the parameter e\, = ei + (a,P,8 = 1) that 
allows factorization of Eq. l l3~35l with 8 = 1. e\ ^ 0. a G [-20,20] and j3 G [-20,20]. 




Fig. 3.9: Imaginary part for the factorization curve of the parameter 
e \+ — j3,5 = 1) that allows factorization of Eq. J3.35> with 5 = 1. e\ ^ 0. 

a e [-20,20] and j8 e [-20,20]. 



Eqs. (I3.37t and (I3.42t representing particular solutions for the GBHE, are the same 
as those obtained by Wang et al. ll27l . 
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3.6 



Conclusion of the chapter 



In this chapter, we apply the same factorization scheme for more complicated 
second order nonlinear differential equations as in Chapter 2. Exact particular 
solutions have been found for a series of nonlinear differential equations with 
applications in physics and biology: the modified Emden equation, the generalized 
Lienard equation, the Duffing-van der Pol equation, the convective Fisher equation, 
and the generalized Burgers-Huxley equation. Also, we display parametric curves 
along which the differential equations under consideration could be factorized. We 
find that the proposed factorization procedure is easier and more efficient than other 
methods used to find particular solutions of second order differential equations. 
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4 



One-parameter supersymmetry for microtubules 



Abstract. The simple supersymmetric model of Caticha 1341 as used by Rosu 1331 to 
describe the motion of ferrodistortive domain walls in microtubules (MTs), is generalized 
to the case of Mielnik's one-parameter nonrelativistic supersymmetry 11 II . By this means, 
one can introduce Montroll double-well potentials with singularities that move along the 
positive or negative travelling direction depending on the sign of the free parameter of 
Mielnik's method. Possible interpretations of the singularity are microtubule associated 
proteins (motors) or structural discontinuities in the arrangement of the tubulin molecules. 



4.1 Introduction 



Based on well-established results of Collins, Blumen, Currie and Ross |29j regarding 
the dynamics of domain walls in ferrodistortive materials, Tuszyhski and collaborators 
1301 13 II considered MTs to be ferrodistortive and studied kinks of the Montroll type 1321 
as excitations responsible for the energy transfer within this highly interesting biological 
context. 

The Euler-Lagrange dimensionless equation of motion of ferrodistortive domain walls as 
derived in 11291 from a Ginzburg-Landau free energy with driven field and dissipation 
included is of the travelling reaction-diffusion type 

y/'+PV - V 3 + V/+a = , (4.1) 

where the primes are derivatives with respect to a travelling coordinate £, = x — vt , p is a 
friction coefficient and a is related to the driven field 1291 . 

There may be ferrodistortive domain walls that can be identified with the Montroll kink 
solution of Eq. ( 14. li 

M ^= ai+ rr5/3l)' (4 ' 2) 

where /3 = [aq, — (X\)j\f2 and the parameters a.\ and are two nonequal solutions of the 
cubic equation 

{y - a x ){y - a 2 ){y - a 3 ) = y/ 3 - y-a . (4.3) 
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4.2 Caticha's supersymmetric model as applied to MTs 



Rosu has noted that Montroll's kink can be written as a typical tanh kink |33 



( ^ 



(4.4) 



where y = (X\ + (X2 = 1 + "'^ . The latter relationship allows one to use a simple 
construction method of exactly soluble double-well potentials in the Schrodinger equation 
proposed by Caticha 1341 . The scheme is a non-standard application of Witten's 
supersymmetric quantum mechanics J9) having as the essential assumption the idea of 
considering the M kink as the switching function between the two lowest eigenstates of the 
Schrodinger equation with a double-well potential. Thus 



</>! = M<j) 



(4.5) 



where fo^ are solutions of <j) Ql + [eo,i — m(^)]^)o,i(^) = 0, and u(%) is the double-well 
potential to be found. Substituting Eq. J4.5i into the Schrodinger equation for the subscript 
1 and substracting the same equation multiplied by the switching function for the subscript 
0, one obtains 

<pi+R M <Po = 0, (4.6) 



where Rm is given by 



M +sM 
2M' 



(4.7) 



and e = £i — £q is the lowest energy splitting in the double-well Schrodinger equation. 
In addition, notice that Eq. J4.6I is the basic equation introducing the superpotential R in 
Witten's supersymmetric quantum mechanics, i.e., the Riccati solution. For Montroll's 
kink the corresponding Riccati solution reads 



sinh(j3£)+2ycosh 2 ( ^ 



(4.8) 



and the ground-state Schrodinger function is found by means of Eq. (14.61 



0o,m(£) = ^Wcosh^jexp^^^exp^-^ cosh(/3 %) 
- 7 j^- 7 sinh(/^f 



(4.9) 



while 0i is obtained by switching the ground-state wave function by means of M. This 
ground-state wave function is of supersymmetric type 



0o,m(4) = <ft)jw(0)exp 



S 



RM{y)dy 



(4.10) 



where 0o,m(O) is a normalization constant. 
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The Montroll double well potential is determined up to the additive constant £o by the 
'bosonic' Riccati equation 



R M - R M + £ 



P 2 (f-l)e 2 e 
4 4j3 2 2 



So 



, L (4^6 + 2( 7 2 + l)ecosh(j3^ ) - 8/3 2 ) cosh(j3^ ) 
-4y (e + £cosh(j3<i; ) - 2/3 2 ) sinh(/3<i; ) 



(4.11) 



Plots of the asymmetric Montroll potential and ground state wave function are given in 
Figs. 4.1 and 4.2 for a particular set of the parameters. If, as suggested by Caticha, one 
chooses the ground state energy to be 



6o 



4 



C 2 + W { '- f) ' 



(4.12) 



then «Af(4) turns into a travelling, asymmetric Morse double-well potential of depths 
depending on the Montroll parameters j3 and y and the splitting e 



U, 



L,R 
0,m 



J3 2 



1± 



ley 



(4.13) 



where the subscript m stands for Morse and the superscripts L and R for left and right well, 
respectively. The difference in depth, the bias, is A m = Uq — Uq = ley, while the location 
of the potential minima on the travelling axis is at 



:L,R 



(2/3) 2 ±2ey 



e(7Tl] 



(4.14) 



that shows that y ^ ± 1 . 



4.3 The Mielnik extension 



We now discuss shortly in this context the Mielnik extension of these results II ll l3l. The 
point is that Rm in Eq. J4.8I > is only the particular solution of Eq. ( 14.1 1> . The general 
solution is a one-parameter function of the form 



d_ 



ln(/ M (|)+A) 



and the corresponding one-parameter Montroll potential is given by 



(4.15) 



(4.16) 



In these formulas, 7m(4) = / ^qm(^)^^ an< ^ ^ ^ s an integration constant that is used 
as a deforming parameter of the potential and is related to the irregular zero mode. The 
one-parameter Darboux-deformed ground state wave function can be shown to be 



0o,M(§;A) = VA(I+r 



(4.17) 
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where \JX{X + 1) is the normalization factor implying that A ^ [0,-1]. Moreover, the 
Mielnik parametric potentials and wave functions display singularities at X s = —Im{<^s)- 
Plots of u M (^;X) and 0o.m(|;A) for A = 1, 10, 100 are presented in Figs. 4.3-4.12 and are 
useful to see the behavior of the singularities and the deformation effect of the X parameter. 
For large values of ±A the singularity moves towards =F°° and the potential and ground state 
wave function recover the shapes of the non-parametric potential and wave function as can 
be seen in Figs. 4.11 and 4.12, respectively. The one-parameter Morse case corresponds 
formally to the change of subscript M — ► m in Eqs. J4. 15i and J4.16> . For the single well 
Morse potential the one-parameter procedure has been studied by Filho 1351 and Bentaiba 
etal 1551 . 

Mielnik's approach leads to singularities in the double- well potential and the corresponding 
wave functions. If the parameter A is positive the singularity is to be found on the 
negative % axis, while for negative A it is on the positive side. Potentials and wave 
functions with singularities are not so strange as it seems (37 1 and could be quite relevant 
even in nanotechnology 1381 where quantum singular interactions of the contact type 
are appropriate for describing nanoscale quantum devices. We interpret the singularity 
as representing the effect of an impurity moving along the MT in one direction or the 
other depending on the sign of the parameter A. The impurity may represent a protein 
attached to the MT or a structural discontinuity in the arrangement of the tubulin molecules. 
This interpretation of impurities has been given by Trpisova and Tuszyhski in non- 
supersymmetric models of nonlinear MT excitations 1391 . 



4.4 Conclusion of the chapter 



In conclusion, the supersymmetric approaches allow for a number of interesting exact 
results in this biological framework and point to a direct connection between Schrodinger 
double-well potentials and nonlinear kinks encountered in nonequilibrium chemical 
processes. MTs are an important application but the procedures described here can be used 
in many other applications. Moreover, the supersymmetric constructions can be used as a 
background for clarifying further details of the exact models. Although it is not so clear 
why one should take a certain type of kink as switching function between the Schrodinger 
split modes, it is interesting that proceeding in this way one will be led to some familiar 
double-well potentials in chemical physics. 
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Fig. 4. 1 : The Montroll asymmetric double-well potential (MDWP) calculated using 
Eq. ETT} for £ = 0. In all figures a, = 1, a 2 = — 1.5, j8 = -2.5/V2, y = -0.5, £ = 0.1. 




Fig. 4.3: The one-parameter Darboux modified MDWP for A = 1. 
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Fig. 4.4: The low-scale left hand side of the singularity. 





Fig. 4. 10: Plot of the integral Im(£ ) that produces the deformation of the potential and 

wave functions. 
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5 



Supersymmetric method with Dirac parameters 



Abstract. In this chapter we first describe a "supersymmetric" one-dimensional matrix 
procedure similar to relationships of the same type between Dirac and Schrodinger 
equations in particle physics that we apply to two problems in classical mechanics and 
quantum mechanics, respectively. In the first case, we obtain a class of parametric 
oscillation modes that we call K-modes with damping and absorption that are connected 
to the classical harmonic oscillator modes through this supersymmetric procedure that is 
characterized by coupling parameters. When a single coupling parameter, denoted by K, 
is used, it characterizes both the damping and the dissipative features of these modes. 
Generalizations to several K parameters are also possible and lead to analytical results. If 
the problem is passed to the physical optics (and/or acoustics) context by switching from 
the oscillator equation to the corresponding Helmholtz equation, one may hope to detect 
the K-modes as waveguide modes of specially designed waveguides and/or cavities. In 
the second case, the same method is presented in a style appropriate for truly quantum 
mechanical problems and an application to the Morse potential is performed. We obtain 
the corresponding nonhermitic Morse problem with possible applications to the diffraction 
on optical lattices. 



5.1 Introduction 



Factorizations of differential operators describing simple mechanical motion have been 
only occasionally used in the past, although in quantum mechanics the procedure led 
to a vast literature under the name of supersymmetric quantum mechanics initiated by 
a paper of Witten |9J. However, as shown by Rosu and Reyes 1401 . for the damped 
Newtonian free oscillator the factorization method could generate interesting results even 
in an area settled more than three centuries ago. In this chapter, we apply some of the 
supersymmetric schemes to the basic classical harmonic oscillator. In particular, we show 
how a known connection in particle physics between Dirac and Schrodinger equations 
could lead in the case of harmonic motion to chirped (i.e., time-dependent) frequency 
oscillator equations whose solutions are a class of oscillatory modes depending on one 
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more parameter, denoted in the following by K, besides the natural circular frequency (Oq. 
The parameter K characterizes both the damping and the losses of these "supersymmetric" 
partner modes. Moreover, we do not limit this study to one K parameter extending it to 
several such parameters still getting analytic results. Guided by mathematical equivalence, 
possible applications in several areas of physics are identified. Moreover, in the final part 
of the chapter the same supersymmetric scheme is used in the context of exactly solvable 
quantum problem of the Morse potential. A nonhermitic version of the Morse problem is 
introduced in this way. 



5.2 Classical harmonic oscillator: The Riccati approach 



The harmonic oscillator can be described by one of the simplest Riccati equation 

u'+u 2 +fcco 2 = 0, K = ±l, (5.1) 

where the plus sign is for the normal case whereas the minus sign is for the up side down 
case. Indeed, employing u = — one gets the harmonic oscillator differential equation 

w" + k-cOqW = , (5.2) 

with the solutions 

W + cos ( coot + <p+) if K = 1 



w b : 



W_sinh(c0ot + <p_) if K" = — 1 



where W± and q>± are amplitude and phase parameters, respectively, which can be ignored 
in the following. 

The particular Riccati solution of Eq. i5.l\ are 

— COt)tan(caot) if K = 1 
fflDCoth(fflot) if K = — 1 . 

It is well known that the particular Riccati solutions enter as nonoperatorial part in the 
common factorizations of the second-order linear differential equations that are directly 
related to the Darboux isospectral transformations 1131 . 
Thus, for Eq. (15. 2> one gets (D t = 4) 

(D, + Up) (D t - Up) w = w" + (-u p - u 2 )w = . (5.3) 

To fix the ideas, we shall use the terminology of Witten's supersymmetric quantum 
mechanics and call Eq. ( 15. 3> the bosonic equation. We stress here that the supersymmetric 
terminology is used only for convenience and should not be taken literally. Thus, the 
supersymmetric partner (or fermionic) equation of Eq. ( 15. 3> is obtained by reversing the 
factorization brackets 

(D, - Up) (D, + Up) w f = w" + (u p - u 2 )w = w" + co f 2 (t)w = , (5.4) 

which is related to the fermionic Riccati equation 

u'-u 2 -co t 2 (t) =0 , (5.5) 
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where the free term (of is the following function of time 

<% 2 (t)=Up-u^ = 
The solutions (fermionic zero modes) of Eq. i5Al are given by 



2 _ J (Oq{— 1 - 2tan 2 ft)()t) if K- = 1 
"p~ U p I tog(l -2coth 2 coot) ifJC=-l 



Wf : 



cos((«Qt) 

tap 
sinh(fl)ot) 



if jc=1 
if JC=-1 



and thus present strong periodic singularities in the first case and just one singularity at the 
origin in the second case. These 'partner' oscillators, as well as those to be discussed in 
the following, are parametric oscillators, i.e., of time-dependent frequency. Moreover, their 
frequencies can become infinite (periodically). In general, signals of this type are known 
as chirps. "Infinite" chirps could be produced, in principle, in very special astrophysical 
circumstances, e.g., close to black hole horizons 1411 . 



5.3 Matrix formulation 



Using the Pauli matrices ay 
equation 







and a x 



1 



i J V 1 

# W=[a y D t + cj x (iup)]W=0, 



, we write the matrix 



(5.6) 



[ ^ ) is a two component spinor. Eq. i5.6\ is equivalent to the following 



decoupled equations 



(iD t +iu p )wi = 
(— iD t + iu p )w2 = . 



(5.7) 
(5.8) 



Solving these equations one gets wi « coq/ cos(ft)of) and W2 « (Oocos(a>ot) for the K = 1 
case and wi <Bo/sinh(<»of) and W2 °= O)osinh(coof) for the K = — 1 case. Thus, we obtain 



W = 



W] 



Wf 

w b 



(5.9) 



This shows that the matrix equation is equivalent to the two second-order linear differential 
equations of bosonic and fermionic type, Eq. j5.2t and Eq. ( I5.4> , respectively, a result quite 
well known in particle physics. Indeed, a comparison with the true Dirac equation with a 
Lorentz scalar potential S(x) 



[-ic7 y D x + a x (m + S(x))]W = EW , 



(5.10) 



shows that Eq. (15. 6> corresponds to a Dirac spinor of 'zero mass' and 'zero energy' in 
an imaginary scalar 'potential' iu p (t). We remind that a detailed discussion of the Dirac 
equation in the supersymmetric approach has been provided by Cooper et al 1142 1 in 1988. 
They showed that the Dirac equation with a Lorentz scalar potential is associated with a 
susy pair of Schrodinger Hamiltonians. This result has been used later by many authors in 
the particle physics context |43 1. 
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5.4 Extension through parameter K 



We now come to the main issue of this work. Consider the slightly more general Dirac-like 
equation 

# K W = [(J y D t + (7 x (iUp + K)]W = KW, (5.11) 

where K is a (not necessarily positive) real constant. On the left hand side of the equation, 
K stands as an (imaginary) mass parameter of the Dirac spinor, whereas on the right hand 
side it corresponds to the energy parameter. Thus, we have an equation equivalent to a 
Dirac equation for a spinor of mass iK at the fixed energy E = iK. This equation can be 
written as the following system of coupled equations 



iD t wi + (iu p + K)wi = KW2 
-iD t W2 + (iu p + K)w2 = Kwi . 



(5.12) 
(5.13) 



The decoupling can be achieved by applying the operator in Eq. ( 15 . 1 21 to Eq. J5 . 1 31 . For 
the fermionic spinor component one gets 



O) 2 



2K 

(1 +2tan 2 <»ot)+i — tanfflot 



Oh 



Wj = for K = 1 



2K 



D 2 Wj +a>o (1 -2coth 2 <Bot) + i — cothcoot 



C0(, 



whereas the bosonic component fulfills 

D t 2 w+ + co ( 2 



, .2K 
1 — l — tancoot 
(Oo 



for K ■ 



for K = 1 



(5.14) 
(5.15) 

(5.16) 



D t 2 w 2 - a>l 



, - 2K u 
1 — i — cothcoot 



w 2 = for k = - 1 



(5.17) 



The solutions of the bosonic equations are expressed in terms of the Gauss hypergeometric 
functions 2F1 



w+(t;a+,j3 + ) = a+zf-^z^aFi 



p + q,p + q-l,2p;-izi 



^ +e -2ip^(p-i) z -(p-^) z (q-l) 2Fl 



q-p,q-p + l,2-2p;--zi 



(5.18) 



and 



r + s,r + s + 1, 1 +2r;-z 3 



s — r + 1 , s — r, 1 — 2r; -Z3 



w 2 (t;a-,j3_) = a-z^z^Fi 
+ j 8_4 1 Z3 r z| 2 Fi 

where the variables z, (/ = 1, ...,4) are given in the following form: 
zj,2 = itan(toot) =p 1, Z3 4 = coth(o)ot) ± 1, 



(5.19) 
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respectively. The parameters are the following: 
1 



1 -■ 



2K 



2K 
(Oo 



1 L .2K 

2 V tOO 



1 -i- 



.2K 
COo 



The fermionic zero modes can be obtained as the inverse of the bosonic ones. Thus 

1 _ 1 



w 



w+(t;a+,j3+) 



w 



w 2 (t;a_,j3_) 



(5.20) 



A comparison of w+ with the common 1/ cost fermionic mode is displayed in Figs. 5.3 
and 5.4. 

In the small K regime, K <C (Oq, one gets 



w+(t;a + ,j3+) » a+z*? 2) zJ 2) 2 F, 2,1,2- £;-^ Zl (t) 



K 1 



K 



Oo 2 
K K 1 



and 



w 2 (t;a_,/3_) wa_z r 3 z| 2 Fi 1,2,2 + i— ;-z 3 (t) 



,1 + — , — ;--zi(t) 
(Oq (Oo (Oq 2 



K 1 



(5.21) 



(Oo 2 



, • K ■ K ■ K 1 M 

1 -l — ,-i — ,-i — ;-z3(t) 

(Oo (Oo (Oo 2 



(5.22) 



Examining the bosonic equations, one can immediately see that the resonant frequencies 
acquired resistive time-dependent losses whose relative strength is given by the parameter 
K. The fermionic equations having time-dependent real parts of the frequency can 
be interpreted as parametric oscillators which are also affected by losses through the 
imaginary part. 



5.5 More K parameters 



A more general case in this scheme is to consider the following matrix Dirac-like equation 
D t 



-i 

1 



1 W iUp + Ki 

1 J\ iu p + K 2 



(o' 








k ) 


U) 



The system of coupled first-order differential equations will be now 

— iD t + iu p + K2 w 2 =KjWi 
iDt + iUp + Kj wi=K 2 w 2 



(5.23) 

(5.24) 
(5.25) 
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and the equivalent second-order differential equations 



D t 2 w ; - 



iAK 



D t w, 



±Dtu p +i(Ki+K 2 )up + (KiK 2 -K 1 K 2 



w,=0, (5.26) 



where the subindex i = 1,2 and AK = Ki — K2. Under the gauge transformation 

1 



w, =Z, exp ( -- 



— iAK 



dr =Z,-(t)e7 1 



{AK 



(5.27) 



one gets 



D t 2 Z, + e,(t)Z,=0, 



where the 'potentials' have the form 

Qi (t) = [ ± D t u p + i (Ki + K 2 )u p + (Ki K 2 - Kj K, ) - u p 2 



1 






-iAK 


~4 





(5.28) 



(5.29) 



Qi,2 are functions that differ from the nonoperatorial parts in Eqs. (I5.55i - d5.33i only by 
constant terms. Indeed, one can obtain easily the following equations. 
For the fermionic spinor component one gets 



D 2 Z+-ft) 2 



, „ , ( AK ) K 1 K 2 -K,K 2 .K!+K 2 
1 + 2tan 2 C0Qt-- i 9 2 +1— -lancet 



OX) 



Z+=0 (5.30) 



for K = 1 , and 



D^ + ft) 2 



, „ i2 (AK) 2 KiK 2 -K 1 K 2 .K!+K 2 , 
1 - 2coth 2 (0 t + £- + 9 1 2 +1— -coth coot 



4co 2 



C0() 



z- = 



for K - = - 1 . 

The bosonic component fulfills 



(5.31) 



D^Zj + col 



1 | (AK) 2 | K^-KX . K!+K 2 



t0() 



tancoot 



Z+=0 



(5.32) 



for fc = 1, and 



D t 2 Z 2 - -coi 



(AK) 2 K t K 2 — K t K 2 . K!+K; 



4« 2 



cothcoot 



Z 2 =0 



(5.33) 



for K" = — 1 . When Ki = K 2 = K one gets the particular case studied in full detail above. 
The more general bosonic modes have the form: 



Zj (t; a+, j3+) = a+ [tan(coot) - i] [tan(fflbt) + i] ^ 
2Fl ^T'^r + 1 ' 1 + 2^ ; 2 (tan ^-) 

£2[ o. t a 2 

+j5+(-iy 7 ^ [tan(coot) - i]"^ [tan(fflbt) + i] 3 ^ 



X 2 F! 



, h 1 , 1 ; -(tanifflnt) — 1 



(5.34) 
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and 



x 2 Fi 



n 3 si4 
Zj (t;a_,j3_) = a_[coth(coot) - l] 5? *[coth(fi) t) + l}^ 

-^r +1 '^T' 1+ 2i ; -2( coth ^- 1 ) 



a 3 n 3 



+j3_(-l) 3S ?4 2S ?[coth(co t)-l] ^[coth(a)ot) + l] : 

X2Fi ^■^ +1 - 1 -dH (c ^ (flW) - 1) 



(ij,. 



(5.35) 



where 



n 3 



Q, = (4co 2 + (K 1 +K 2 ) 2 +4[(K 1 +K 2 )a^-K 1 K 2 ] 
H 2 = (4co 2 + (Ki +K 2 ) 2 -4[(K! +K 2 )o) + K' 1 K 2 ] 



4ojg - (Ki +K 2 ) 2 -4[i(Ki +K 2 )<bd- ^K 2 
4ojg - (Ki + K 2 ) 2 + 4[i(Ki + K 2 )ft)o + KjKj 



1/2 



1/2 



1/2 



5.6 Possible applications of the K-modes 



5.6.1 Waveguides. 

In view of the correspondence between mechanics and optics, one can also provide 
an interpretation in terms of the Helmholtz optics for light propagation in waveguides 
of special profiles. The supersymmetry of the Helmholtz equation has been studied 
by Wolf and collaborators 1441 . To get the waveguide application, one should switch 
from the temporal independent variable to a spatial variable t — > x along which we 
consider the inhomogeneity of the fiber whereas the propagation of beams is along 
another supplementary spatial coordinate z. Thus, we turn the equations d5.55t - d5~.56t 
into Helmholtz waveguide equations of the type (we take c = 1) 

[d 2 + d 2 + O) 2 n 2 (x)]<p(x,z) =0 , (5.36) 

where the modes (p(x,z) can be written in the form w 1 2 (x)e~ lkzZ for a fixed wavenumber k z 
in the propagating coordinate that is common to both wave functions and the index profiles 
correspond to two pairs of bosonic-fermionic waveguides and are given by 

2K 2K 

n£(x)~ 1 -i — tan(kox) , n 2 (x) (1 + 2tan 2 ((0ox)) - i— tan(k x) , (5.37) 

ko k 

and 

nl(x) 1 -i^coth(kox) , n 2 (x) ~ 1 -2coth 2 (k x) +i^coth(k x) , (5.38) 

ko k 

respectively. In our units ko = COq. Eqs. 15.371 . 15.381 can be obtained from Riccati 
equations of the type (c^ 1) 

fflgn? b (x)/c 2 = k 2 =F R x - R 2 , (5.39) 
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where R(x) are Riccati solutions directly related to the Riccati solutions discussed in the 
previous sections. 

According to Chumakov and Wolf 11441 a second waveguide interpretation is possible 
describing two different Gaussian beams, bosonic and fermionic, whose small difference in 
frequency is given in terms of a small parameter e (wavelength/beam width), propagating 
in the same waveguide. In this interpretation, the index profile is the same for both beams. 
For illustration, let us take the normal oscillator Riccati solution in the space variable 
x, i.e., tankox that we approximate to first order linear Taylor term kox. Then, the two 
beam interpretation leads to the following Riccati equation (for details, see the paper of 
Chumakov and Wolf) 



An almost exact, up to nonlinear corrections of order e 2 and higher, supersymmetric 
pairing of the z wavenumbers (propagating constants) occurs, except for the 'ground state' 
one. As noted by Chumakov and Wolf, supersymmetry connects in this case light beams 
of different frequencies but having the same wavelength in the propagation direction z. 
This approach is valid only in the paraxial approximation. Therefore, one should know the 
small x behavior of the K-modes in order to hope to detect them through stable interference 
patterns along the waveguide axis. 

5.6.2 Cavity physics. 

Another very interesting application of the K-modes in a radial variable could be 
Schumann's resonances, i.e., the resonant frequencies of the spherical cavity provided by 
the Earth's surface and the ionosphere plasma layer 1451 . The Schumann problem can be 
approached as a spherical Helmholtz equation [V 2 + k 2 ](j) = with Robin type (mixed) 
boundary condition ^-\s = C(ffl)0s, where C(ffl) is expressed in terms of the skin depth 
8 = \j2/{ii c O(0) of the conducting wall, fX c is its permeability and a is its conductivity. 
The eigenfrequencies fulfilling such boundary conditions can be written as follows 



where / is a complicated expression in terms of skin depths and surface and volume 
integrals of Helmholtz solutions with Neumann boundary conditions ^ \$ = 0. It is worth 
noting the similarity between these improved values of Schumann's eigenfrequencies and 
the K-eigenfrequencies. Moreover, using the Q parameter of the cavity, one can write 
Eq. J5.4H in the form 



This form shows that the modification of the real part of co leads to a downward shift of 
the resonant frequencies, while the contribution to the imaginary component changes the 
rate of decay of the modes. 

We point out that Jackson mentions in his textbook that the near equality of the real and 
imaginary parts of the change in co 2 is a consequence of the employed boundary condition, 
which is appropriate for relatively good conductors. Thus, by changing the form of C(co) 
that could result from different surface impedances, the relative magnitude of the real 
and imaginary parts of the change in co 2 can be made different. It is this latter case that 
corresponds better to the K-modes. 



0) 2 , 2 n 2 (x) - co 2 n 2 (0) = T k -kgx 2 (l =p e) . 



(5.40) 




(5.41) 




(5.42) 
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5.6.3 Crystal models. 

There is also a strong mathematical similarity between the K-modes and the solutions 
of Scarf's crystal model [46:1 based on the singular potential V(x) = — Vocosec 2 (7Tx/a), 
where a is an arbitrary lattice parameter. For this model the one-dimensional Schrodinger 
equation has the form 



where f(x) = sm(nx/a) corresponds to the z;(f) functions, and s and A corresponding to 
-p and -q, respectively, are related to the potential amplitude and energy spectral parameter. 
Thus, by turning the K-oscillator equations into corresponding Schrodinger equations, one 
could introduce another analytical crystal model with possible applications in photonics 
crystals. 

5.6.4 Cosmology. 

Rosu and Lopez-Sandoval apply the K-mode approach to barotropic FRW cosmologies 
1471 . K- Hubble cosmological parameters have been introduced and expressed as 
logarithmic derivatives of the K-modes with respect to the conformal time. For K — > 
the ordinary solutions of the common FRW barotropic fluids have been obtained. 
It is also worth noticing the analogy of the nonzero K oscillator case with the phenomenon 
of diffraction of atomic waves in imaginary crystals of light (crossed laser beams) |48 1. In 
fact, the K parameter is a counterpart of the modulation parameter Q introduced by Berry 
and O'Dell in their study of imaginary optical gratings. Roughly speaking, the nonzero 
K modes could occur in an imaginary crystal of time that could occur in some exotic 
astrophysical conditions. 




(5.43) 



For < x < a/ 2, the general solution is 



¥ =[f(x)]^ 2 F l [I + I(, + A),i + i(,-A);l+,;/ 2 w] + 



[/(*)]KFi [I-I( J -A),i-i(, + A);l-,;/ 2 (x) 



(5.44) 
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Fig. 5.1: The real part of the bosonic mode (y; j, j) for t £ [0, 10] and K £ [0,4] 




Fig. 5.2: The imaginary part of the bosonic mode wj(y; j, 5) for t £ [0, 10] and 

KG [0,4]. 




Fig. 5.3: The real part of the bosonic mode wj(y; A, A) for t £ [0,20] and K = 0.01. 




Fig. 5.4: The imaginary part of the bosonic mode wj(y; I, i) for t £ [0,20] and 

K = 0.01. 
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Fig. 5.5: The real part of the bosonic mode w^(y; j, j) for t € [0,20] and K = 2. 
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Fig. 5.6: The real part of the bosonic mode w^(y; |, ^) for t G [0,20] and K = 2 in the 

vertical strip [-0.5, 0.5]. 
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Fig. 5.7: The imaginary part of the bosonic mode (y; \, k) f° r t € [0,20] and K = 2. 
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Fig. 5.8: The fermionic zero mode — 1/ cost, (red curve), and the real part of — l/w^, 

(blue curve), for K = 0.01. 




Fig. 5.9: The fermionic zero mode — 1/cost, (red curve), and the imaginary part of 
— 1 / w 2, (blue curve), for K = 2. 



5.7 Quantum mechanics with Riccati nonhermiticity 



We have elaborated in the previous sections on an interesting way of introducing imaginary 
parts (nonhermiticities) in second order differential equations starting from a Dirac - 
like matrix equation 1471 1491 . The procedure is a complex extension of the known 
supersymmetric connection between the Dirac matrix equation and the Schrtidinger 
equation. A detailed discussion of the Dirac equation in the supersymmetric approach 
has been provided by Cooper et al. 1421 1431 in 1988, who showed that the Dirac equation 
with a Lorentz scalar potential is associated with a susy pair of Schrodinger Hamiltonians. 
In the supersymmetric approach one uses the fact that the Dirac potential, that we denote 
by S, is the solution of a Riccati equation with the free term related to the potential function 
U in the second order linear differential equations of the Schrodinger type. 
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Indeed, writing the one-dimensional Dirac equation in the form 



[ap + j5m + pR(x)]y(x) = Ey{x) (5.45) 
where c = Ti = 1, p = —id/dx, m (> 0) is the fermion mass, and R(x) is a Lorentz scalar. 



The wave function y is a two-component spinor ( ) and the Pauli matrices a and /3 
are the following 



i o 1 ) and ff * = (l 

Writing the matrix Dirac equation in coupled system form leads to 

[D x + m + R]y 1 =Eyf 2 (5 .46) 

[-D x + m+R]y/2=Ey i (5.47) 
By decoupling one gets two Schrodinger equations for each spinor component, respectively 

HiXj/i = [-D 2 x + Ui\ Wi = £Wi , £ = E 2 - m 2 , (5.48) 

where i = 1,2, and 

Ui (x) = (R 2 + 2mR =F dR /dx) . 
One can also write factorizing operators for Eqs. J5.48> 

A ± = ±D x + m+R (5.49) 

such that 

ill ill 

H l= A-A + ~-, H 2 =A+A--- (5.50) 

However, we have employed the method for the case of the classical harmonic oscillator, 
which is the very specific situation in which the Dirac mass parameter that we denoted 
by K was treated as a free parameter equal to the Dirac eigenvalue parameter E. This 
is equivalent to Schrodinger equations at zero energy, e = 0. On the other hand, it is 
interesting to see how the method works for negative energies, i.e., for a bound spectrum 
in quantum mechanics. Here we briefly describe the method and next apply it to the case 
of Morse potential. 



5.8 1 Complex extension with a single K parameter 

We consider the slightly different Dirac-like equation with respect to Eq. J5.45> 

# K W= [cr y D x + cT x (iR + K)]W = KW , (5.51) 

where K is a (not necessarily positive) real constant. In the left hand side of the equation, 
K stands as a mass parameter of the Dirac spinor, whereas on the right hand side it 
corresponds to the energy parameter. R is an arbitrary solution of the Riccati equation 
of the Witten type |9| 

R'±R 2 = u, (5.52) 
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where u is the real part of the nonhermitic potential in the Schrodinger equations we get. 
Thus, we have an equation equivalent to a Dirac equation for a spinor W = 

(Wf \ 
w J of mass K at the fixed energy E = K but in a purely imaginary potential (optical 

lattices). This equation can be written as the following system of coupled equations 

iD x ^i + (iR + K)^i=K^ ) (5.53) 



-iD x </) o + (iR + K)</. o = K0 1 



(5.54) 



The decoupling of these two equations can be achieved by applying the operator in 
Eq. J5.54I to Eq. d5 .53t . For the fermionic spinor component one gets 



R D X R i2KR 



whereas the bosonic component fulfills 



D 



R 2 + D x R-i2KR 



01=0 



00 = 



(5.55) 



(5.56) 



This is a very simple mathematical scheme for introducing a special type of nonhermiticity 
directly proportional to the Riccati solution. 



5.9 Complex extension with parameters K and K'. 



A more general case in this scheme is to consider the following matrix Dirac-like equation 
D x 



-i 

1 



1 

1 



iR + K 
iR + K 



W] 

Wi 



K 






V o 







The system of coupled first-order differential equations will be now 



iD x + iR + K 



W2 = K W] 



iD x + iR + K wi =Kw 2 
and the equivalent second-order differential equations 

D x 2 w;+ [±D x R + 2iKR+(K 2 -K 2 )-R 2 



W; = , 



(5.57) 

(5.58) 
(5.59) 

(5.60) 



where the subindex i = 1,2 refers to the fermionic and bosonic components, respectively. 
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5.10 Application to the Morse potential 



This potential is frequently used in molecular physics in connection with the disassociation 
of diatomic molecules. In this case, the Riccati solution is of the type 

R=A-Be-" X , (5.61) 

Therefore, the second-order fermionic differential equation will be 

D 2 x wi + [-(Be- 2ax -Cie-" x ) +{K 2 -K' 2 )-A 2 

+ 2iK(A-B<T ax )] wi=0 (5.62) 

where B = B 2 , and C\ = B(2A + a). 

The solution is expressed as a superposition of Whittaker functions 

wi = a^ x l 2 M KupL (^e"^ +j3ie M / 2 W^, M ( 5 - 63 ) 

^ = i(2 + f-/f)and Al = §(4#-/f) 1/2 . 
The bosonic equation reads 

D 2 w 2 + [- (Be- 2ax - C 2 ^ ax ) + (K 2 - K' 2 ) -A 2 

+ 2iK(A - Be-" x )] w 2 = (5.64) 

where B = B 2 , and C 2 = B(2A - a). 

The solution is a superposition of the following Whittaker functions 

w 2 = a 2 e" x / 2 M K2 ^ (™ e-^J + fe^W^ e"^ (5.65) 
where K 2 = ^ (2 — f — i^j-) and the ji subindex is unchanged. 

If we now place ourselves within the quantum mechanical (hermitic) Morse problem we 
should take j3 2 = and K = in order to achieve the exact correspondence with the bound 
spectrum problem and eliminate the nonhermiticity. Moreover, the following well-known 
connection with the associated Laguerre polynomials 

M^ n+ ^(y)=y^-y/ 2 LP(y), y = ™ e~« (5.66) 

can be used in our case with the following identifications 

P K' , , 

? = — . K' = {A-an) 

2 a 

Then we can write the solution of the hermitic bosonic problem in the well-known form 



w 2 ,n{y) = a 2 ( -V y-«-"e-y/ 2 L 2 n { °- n) (y) . (5.67) 
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If we want to approach the nonhermitic problem we define by analogy with Eq. J5.66i 

M^{y)=y^-y/ 2 L 2 l , (y) , (5.68) 

where K% and jj, are the complex parameters mentioned before and the symbol 
corresponding to the associated Laguerre polynomial representing now a Laguerre-like 
function introduced by definition through Eq. (I5.68> . The wave function of the nonhermitic 
problem can be written as follows 

W2 <nonhem iy) = «2 (—) 1 /e-W" H , (y) • (5.69) 

For the nonhermitic fermionic problem, the formulas are similar with the replacement of 
K 2 by Kj . 




Fig. 5.10: Real part of the bosonic wave function w 2 in the range x € [0, 3] and 

Ke [0,2]. 
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Fig. 5.11: Imaginary part of the bosonic wave function W2 in the range x £ [0, 3] and 

Ke[0,2]. 




Fig. 5.12: Real part of the fermionic wave function w\ in the range x <G [0, 3] and 

Ke[0,2]. 




Fig. 5.13: Imaginary part of the fermionic wave function w\ in the range x E [0, 3] and 

Ke[0,2]. 
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5.11 Conclusion of the chapter 



By a procedure involving the factorization connection between the Dirac-like equations 
and the simple second-order linear differential equations of harmonic oscillator type, 
a class of classical modes with a Dirac-like parameter describing their damping and 
absorption (dissipation) has been introduced in this chapter. While for zero values of the 
Dirac parameters the highly singular fermionic modes are decoupled from their normal 
bosonic harmonic modes, at nonzero values a coupling between the two types of modes is 
introduced at the level of the matrix equation. These interesting modes are given by the 
solutions of the Eqs. (I5.55i - J5".33i and in a more general way by Eqs. (I5.27i . (5.34)-(5.35) 
and are expressed in terms of hypergeometric functions. Several possible applications 
in different fields of physics are mentioned as well. Finally, similar to the fact that the 
PT quantum mechanics can be considered as a complex extension of standard quantum 
mechanics, we notice that what we have done here is a particular type of complex extension 
of the classical harmonic oscillator. The complex supersymmetric extension introduced in 
the first sections has been applied to exactly solvable quantum Morse problem. The bosonic 
and fermionic wave functions have been obtained in explicit form. This complex extension 
could have applications to the diffraction of diatomic molecules on optical lattices (systems 
of laser beams). 
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Preliminary remarks on Part II 



The following remarks are pointed out in order to describe the relationship between Part 
I and Part II of the thesis work, where factorization methods for nonlinear ODE and 
synchronization of a neuronal ensemble through feedback methods, respectively, have 
been developed. The study of many biological systems from a mathematical point of 
view is very stimulating because of the possibilities to forecast the dynamic behavior of 
such biological systems. An important example is the neuronal dynamics that governs 
many living organisms. The neuronal dynamics performs the processing of biological 
information by means of transmitted signals. The pulse propagation along a nerve axon 
as described by the FitzHugh-Nagumo equation that was solved in Section 2.4 through 
factorization methods is given by a travelling signal of the kink type. The kink solution 
represents the transition between two stable equilibrium states, and the "level change" 
travelling waves are transmitted along the neuron axon 1501 . It is a very interesting 
challenge to study the synchronization dynamics for a minimal ensemble of two neurons. 
However, the kink type solutions that have been considered by us could not describe 
the most realistic dynamical behavior for the neuronal ensemble. In fact, the FitzHugh- 
Nagumo neuron model can be shown to be an approximation (see, e.g. |50|) of the 
widely known Hodgkin-Huxley (HH) neuron model |69], Therefore, the problem of 
neurons synchronization has been addressed in the more realistic case of the HH neuron 
model. Although theoretical attempts have been made to obtain analytical solutions of 
the travelling wave type for the HH systems (see, e.g. 1511 1521 ). it is an easier task to 
achieve numerical results to study the synchronization dynamics. That is why we do not 
pay attention to the HH travelling waves in the second part. The nonlinear control theory 
is a very suitable way to study and search for the intrinsic mechanisms underlying the 
synchronization phenomena. In the next Chapters 7 and 8, the problem of synchronizing a 
minimal ensemble of two HH neurons is stated, and results on its synchronized dynamics 
are achieved by implementing adequate feedback schemes. 
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Synchronization of chaotic dynamics and neuronal 

systems 



Abstract. In this chapter, the synchronization phenomena and the concept of chaos, their 
importance in natural process and engineering systems are reviewed. A brief overview of 
synchronization methods for the control of chaos and its applications in biological systems 
is presented. Neuronal synchronization activity and its role on brain dynamics is also 
discussed in order to state the problem of neuronal synchronization employing nonlinear 
control theory tools. In addition, the dynamical model for the Hodgkin-Huxley (HH) 
neurons is described. 



7.1 Introduction 



Synchronization phenomena are very important processes occurring in nature and often 
produced as a desired behavior in engineering systems. In very general terms, the 
synchronization of coupled systems means "to share time or events". It refers to the way in 
which networked elements, due to their dynamics, communicate and exhibit collective 
behavior 1531 . Some examples are the observed synchronized flashing of fireflies, 
synchronization of cells in a beating heart, the quantum synchronization in superfluidity 
and superconductivity, in the phenomena related to Josephson tunneling. Other important 
examples are the generated synchronization in computer chips, communication systems 
and global positioning systems. 

On the other hand, the disquieting question about the exact forecast of the evolution in time 
of diverse systems produced the discovery of the existence of chaos. The concept of chaos 
usually refers to the issue of whether or not it is possible to make long-term predictions 
about the behavior of a system. There exist several mathematical definitions of chaos, 
however, all of them express the property of high sensitivity to the initial conditions. Such 
a characteristic implies that two even arbitrarily close trajectories, separate exponentially 
in the course of time. A deterministic system is said to be chaotic if necessary requirements 
of nonlinearity and dimensionality (of at least three) are characteristic for that system 1 54 1 . 



66 



During the last fifteen years there have been an increased interest in studying chaotic 
systems, since the proved fact that chaos "cannot be forecasted but it can be controlled" 
1 55 1. The issue of the control of chaos is of interest for both theorist and control engineers. 
Interest arises because of many observations show that the chaotic behavior is common in 
nature, for instance, chaotic dynamics can be found in meteorology, plasma physics, heart 
and brain of living organisms; and experimental results dealing with the control of chaotic 
systems lead us to practically an unlimited amount of technical applications in mechanical 
and space engineering, electrical and electronic systems, communication and information 
systems, etc. 1561 . 

Two applications of the control of chaos have been widely studied for the past few year: 
the control and use of chaos for communication systems, and the synchronization (and 
suppression) of chaotic dynamics for several communication schemes B54II57I . 
In the following section a review of synchronization methods concerning the control of 
chaotic dynamics and its applications in biological systems is presented. Then we focus 
on synchronization of a neuronal system comprised of two isolated HH neurons, whose 
dynamical model is also described in final section. In Chapter 8, the obtained results 
employing the mathematical tools provided by the nonlinear control theory 1581 . that can 
be applied in order to show the way the HH neurons can synchronize are presented. We 
find that isolated neurons unidirectionally couple and synchronize through feedback action. 
In addition, robust synchronization dynamics is obtained by implementing a dynamic 
compensator. 



7.2 Synchronization methods for the control of chaos 



Many approaches have been proposed to control chaotic dynamics of systems. Open loop 
or non feedback methods and closed loop or feedback methods 1541 1531 1561 have been 
developed in order to produce the desired behavior in chaotic systems. Two basic problems 
concerning the control of chaotic dynamics through feedback methods are identified: 
synchronization and suppression of chaos. Suppression of chaos consists in stabilization 
of the system around regular orbits or equilibrium points. The chaos synchronization 
problem has the characteristic that the receiver (slave) system must track in some sense 
the trajectories of the sender (master) system 1591 . Because of uncertainties may appear in 
the chaos control problem, adaptive schemes are implemented in order to achieve robust 
synchronization dynamics. Synchronization can be achieved for identical chaotic systems 
with (obviously) different initial conditions [60|, however, researchers have found that 
non identical systems synchronize when adequate developed feedback schemes are applied 

From the standpoint of the geometrical control theory, the synchronization problem can 
be seen as a stabilization problem. For a defined synchronization error x e = x e ,u — x e ,s, 
where x e represents the difference between the master and slave system states, and 
x S R d , e = l,2,...,d, there exists a synchronization error system H6II whose trajectories 
exponentially converge to zero under a feedback control action; consequently, the master 
and slave systems unidirectionally couple. 

In the following chapter, the statement of the synchronization problem and its solution 
using the geometrical control theory for a proposed chaotic neuronal system is explained 
in detail. Also, because of appearing uncertain system states, i.e., non accurately 
measured states, an adaptive scheme is implemented via construction of a state observer or 
uncertainty estimator that guarantees robust synchronization. 
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7.3 Applications of synchronization methods in biological systems 



Synchronization methods for the control of chaos where feedback action is implemented 
have a wide variety of technical and scientific applications. Among the main technical 
applications of chaos synchronization can be found a diversity of schemes for communica- 
tion systems. The scientific applications are directed to study properties, regularities, and 
mechanisms of the behavior of physical, chemical and biological systems. Moreover, very 
interesting results are obtained when the control methods are applied to experimental sys- 
tems. Examples of suppression and synchronization of chaos can be found in biomechani- 
cal systems, medicine, biology and ecology: design of feedback pacemakers, suppression 
of oscillating epileptiform activity in neural networks, control of population dynamics in 
plankton and other biological species, etc. 

The main interest in Part II is to apply feedback synchronization methods to a chaotic 
neuronal system. Two HH neurons are regarded as a system of two dynamical subsystems, 
with the aim to show that synchronized dynamics is achieved through feedback strategies. 
Synchronization and suppression of chaos, using the tools of control theory, could provide 
insights to understand and show their relevance in the processing of biological information 
in neural networks. 



7.4 Synchronized dynamics of neurons 



Synchronization of neuronal activity patterns (action potentials) is a fundamental 
topic in the modern research of brain dynamics. Experimental evidence reveals that 
synchronization phenomena are basic for the processing of biological information. It has 
been demonstrated |62J that large ensembles of neurons whose functionality is related to 
visual perception synchronize their oscillatory activity (in the gamma frequency range, 40- 
60 Hz) when stimulated. Some researchers consider that neuronal synchronization allows 
the brain to solve the so-called binding problem 1631 . Take for example a car: it may be 
characterized by its shape, color, emitted noise, and so on. All these features are processed 
in different parts of the brain, however we conceive the car as a single entity because of 
still unknown binding procedures. More recent studies suggest that synchronization is a 
basic mechanism for consciousness 1641 1651 l66l . Also, Parkinsonian tremor and epileptic 
seizures are closely related with this mechanism 1631 1671 . 

Recently, it was shown that two coupled living neurons synchronize their activity patterns 
when depolarized by an external current 1 68 1. However, the whole underlying mechanisms 
are not completely understood. From a theoretical point of view, neurons are considered 
as nonlinear oscillators. A lot of theoretical studies have been carried out to investigate 
the dynamics of single neurons and neural networks. The most employed and realistic 
neuron models are the Hodgkin-Huxley and the FitzHugh-Nagumo systems 1631 . that have 
been used in detailed studies of neuron behavior under external forcing. Moreover, within 
these models one can consider the change of dynamical parameters and implications on the 
neuronal activity, for instance, the strength of synaptic conductance, and intrinsic or added 
noise. The neuronal synchronization problem is addressed regarding diffusive coupling, or 
modelling unidirectionally coupled master-slave systems. Many people believe that control 
theory could be very useful to address the problem of synchronization in ensembles of 
neurons. The mathematical tools provided by the nonlinear control theory, can suitably be 
applied to the mathematical model of a system comprised of two noiseless HH neurons. 
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7.5 The Hodgkin-Huxley model of the neuron 



The brain is the most complex system known to us. Understanding the way it works and its 
structure has been a very interesting research issue during many decades. The basic units 
which integrate the brain tissue and, in general, any nervous system, are the nerve cells 
or neurons. The transmission of signals and the processing of biological information are 
carried out through complicated interactions between large ensembles of neurons. External 
and internal stimuli generate the biological information which propagates from the sensory 
sites to specific areas of the brain, for instance, the visual, olfactory and auditive perception, 
and the conscious sensory-motor activity. 

Most of the nerve cells generate a series of voltage spiking sequences called the membrane 
action potential, in response to external stimuli performing the information processing. 
These pulses of the action potential originate at the cell body and propagate down the 
axon at constant amplitude and velocity. They can be transmitted via synaptic coupling to 
another nerve cell which is stimulated by the corresponding current. The electric behavior 
of the cell axon membrane is described by the net ion flux through a great amount of 
potassium and sodium ionic channels; each ion passing from the inner (outer) to the outer 
(inner) side of the cell. It is known that potassium and sodium channels are composed 
of four independent gates, which can be in a permissive or non permissive state. The 
potassium ions cross the membrane only through channels that are specific for potassium. 
If the four gates of a potassium channel are in permissive state, then the channel is open 
and potassium ions flow through it. The sodium ions cross the membrane only through 
channels that are specific for sodium. Three of the gates for a sodium channel are activation 
gates, and one is inactivation gate. All of them must be in permissive state to allow sodium 
ions to cross the sodium channel. 

Several neuron models have been proposed to describe the dynamics of the action 
potentials. However, the most widely used is still the realistic HH neuron model 1691 . In 
the early 1950's, Hodgkin and Huxley developed and published a series of investigations 
where they studied the electrophysiology of the squid giant axon. Their results allow them 
to calculate the total membrane current as the sum of the potassium and calcium ionic 
channels currents and the capacitive current, 

r m (t)=Iionic(t)+C m ^}p-. (7.1) 
dt 

In addition, based on the large number of realized experiments they postulated a phe- 
nomenological model that turned itself into a paradigm for the generation of the action 
potential in the squid axon. We follow next the recent discussion in the book of C. Koch 
1701 for a compact presentation of the main statements of the HH model: 

1. The action potential involves two major voltage-dependent ionic conductances, a sodium 
conductance g^ a and a potassium conductance gK- They are independent from each other. 
A third, smaller "leak" conductance g/ does not depend on the membrane potential. The 
total ionic current flowing is given by the following equation 

Iionic(t) = ha +Ik + keak- (7-2) 

2. The individual ionic currents 7,(f ) are linearly related to driving potential via Ohm's law, 

I,(f)=g,{V(t),t)(V(?)-E,) (7.3) 
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where the ionic reversal potential £, is given by Nernst's equation for the appropriate ionic 
species. Depending on the balance between the concentration difference of the ions and 
the electrical field across the membrane separating the intracellular cytoplasm from the 
extracellular milieu, each ionic species has an associated "ionic battery". Conceptually, 
there exist an equivalent electrical circuit to describe the axonal membrane. 
3. Each of the two ionic conductances is expressed as a maximum conductance, g^ a 
and gK, multiplied by a numerical coefficient representing the fraction of the maximum 
conductance actually open. These numbers are functions of one or more Active gating 
particles Hodgkin and Huxley introduced to describe the dynamics of the conductances. 
In their original model, they talked about activating and inactivating gating particles. Each 
gating particle can be in one of two possible states, open or close, depending on time and on 
the membrane potential. In order for the conductance to open, all of these gating particles 
must be open simultaneously. The entire kinetic properties of their model are contained in 
these variables. 

The gating particles, also known as gating variables, Hodgkin and Huxley presented are 
usually denoted by n, m and h. They are the same as the currently known ion channel gates. 
The following set of four coupled nonlinear differential equations represents the complete 
HH neuron dynamical model 1691 : 

C m ^- = Iext-gKn 4 {V-V K )-gNam 3 h(V-V Na )-gi(V-V Na ), (7.4) 



dt 
dn 

dt 
dm 

~dt 
dh 

dt 



an(V){l-n)-p n (V)n, (7.5) 
a m (V)(l-m)-p m (V)m, (7.6) 
a„(V)(l-h)-p h (V)h, (7.7) 



where V represents the membrane potential, n is the probability of any given potassium 
channel gate being in the permissive state (activation of the potassium flow current), m is 
the probability of any given activation sodium channel gate being in the permissive state 
(activation of the sodium flow current), and h is the probability of any given inactivation 
sodium channel gate being in the permissive state (inactivation of the sodium flow current). 
Gating variables are dimensionless and within the range [0,1]. C m is the membrane 
capacitance, gK, gNa and g/ are the maximum ionic and leak conductances, while Vk, Vnu 
and Vi stand for the ionic and leak reversal potentials. The external stimulus current can be 
modelled by the term I ext , usually a tonic or periodic forcing. The respective differential 
equations for n, m and h, describe the transition from open to closed states for the gating 
variables. The explicit form of the functions ttj(V) and fij(V) (j = n,m,h) in Eqs. (17. 51 - 
il.li is given as follows 1691 1701 . 

a "= { ^TiO)Z-iy A^0.125exp(y/80); (7.8) 

a ^ { exp[(y+ V 25 + )/?0]-l} ' A^4exp(y/18); (7.9) 

a /7 =0.07exp(V/20), ft = {ap[(y + M)/10] + 1} ' (7 ' 10) 



Also, nominal values for the system parameters can be found in I?69I I70I . 
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Unidirectional synchronization of Hodgkin-Huxley 
neurons 



Abstract. Synchronization dynamics of two noiseless Hodgkin-Huxley (HH) neurons 
under the action of feedback control is studied. The spiking patterns of the action potentials 
evoked by periodic external modulations attain synchronization states under the feedback 
action. Numerical simulations for the synchronization dynamics of regular-irregular 
desynchronized spiking sequences are displayed. The results are discussed in context of 
generalized synchronization. It is also shown that the HH neurons can be synchronized in 
face of unmeasured states. 



8.1 Introduction 



For several decades many attempts have been addressed to understand the processing 
of biological information in single neurons and neural networks. Experimental reports 
171117211731 suggest that the synchronization plays a very important role in the processing 
of information by large ensembles of neurons. Recently, it has been demonstrated that a 
minimal ensemble of two coupled living neurons fire synchronized spiking activity when 
depolarized by an external DC current 1741 . However, total neural mechanisms underlying 
synchronization are not well understood yet. The Hodgkin-Huxley neurons are usually 
used as realistic models of neuronal systems, for studying neuronal synchronization. 
Some theoretical approaches investigate the synchronization phenomena considering 
diffusive coupling and the influence of intrinsic noise as a promoter of neuronal activity 
1751 1761 . and studying the synchronization dynamics related to the rhythmic oscillations 
phenomena (theta and gamma frequency rhythms) in neurons of localized areas of the 
brain 1771 l78l 1791 . In addition, the forcing of HH neurons by external stimulus has been 
widely studied [80 18 1 1 l82l 1831 l84l for tonic or periodic currents that trigger the action 
potential displaying spike activity and refractory dynamics. 

On the other hand, synchronization of chaotic systems is a relatively recent phenomena 
1 60 1 which can be understood from nonlinear geometrical control theory [59, 61 1. In this 
chapter we study the synchronized behavior of two silent (i.e., the autonomous HH systems 
exhibit fixed point dynamics) 1811 HH neurons, proposing an unidirectionally coupled 
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synchronization system. In the chaos synchronization problem the trajectories of a slave 
system must track, in some sense, the trajectories of a master system even though slave 
and master systems may be different. The obtained results contribute in the theoretical 
framework of neurons synchronization, and relates the phenomena to the well-posed 
concept of generalized synchronization (GS) 1851 . Because of uncertain states cannot be 
accurately measured in practice, they are not available to do control, for instance, the ionic 
channels activation. Then, an approach for robust synchronization via construction of an 
uncertainty estimator, is implemented. 

The organization of this chapter is as follows. In Section 8.2 the HH neuronal systems are 
described. The statement of the problem for unidirectionally coupling synchronization of 
HH neurons is given in Section 8.3. The HH neuronal synchronization dynamics obtained 
through a stabilizing control law is studied in Section 8.4. In Section 8.5 the generalized 
and robust synchronization are discussed, and a final conclusion section is given. 



8.2 The Hodgkin-Huxley system redefined 



We redefine the HH system of equations in order to state the synchronization problem of 
two HH neurons. Let x^m and X{g (i = 1,2,3,4 and subscripts M,S stand for the master 
and slave system, respectively) be the four variables V, n, m and h in each system. As 
the meaning of synchronous behavior is "to share time or events", we shall consider two 
HH neurons modelled by Eqs. \1 A\ - RTn\ . Thus, the master system is represented by the 
following set of equations: 

X 1,M = 1 / C,n M [Iext M ( f ) ~ g^M X 2,M ( x l M ~ V K M ) 

-8Na M x 3,M X4 > M ( Xl > M _ V Na M ) ~gl M ( X l-M _ V l M )1 > ( 8 - *) 

X2.M = «„(xi, M ) (1 -X2.m) - Pn(xi iM )X2,M > ( 8 -2) 

X3,m = (X m {x\,M)( 1 -x 3M )- f5 m (x LM )x 3M , (8.3) 
X4,M = OC h (xiji) (1 - X^m) - Ph(xiJi)X4J{ , (8-4) 

and the slave system is proposed to be governed by the equations: 

*1,S = l/Cms [lext s (t) -8K s X2 : s(xi,S-Vk s ) 

-gNa s xl tS X4,S {xi.S - Vnc, s ) ~gl s (x\ t S ~ Vfc)] + U, (8.5) 
X2.S = Ci„(x LS )(l -X 2> s)- Pn{Xl,s)X2,S, (8-6) 

xi,s = a m (xi t s){l-x3,s)-p m (xi tS )x3,s, (8.7) 
U,s = oc h {x LS )(l -jc 4 , s )-j3/,(x IjS )jc 4iS , (8.8) 

where the added term u in Eq. J8.5I represents a feedback synchronization force. Nominal 
values are considered for parameters of the master system: C mu = 1 fiF /cm 2 , gK M = 36 
mD.- 1 /cm 2 , g NctM = U0mQ.- l /cm 2 , g, M = 0.3 mCT x /cm 2 , V Km = 12 mV, V NaM = -115 
mV and V[ M = —10.613 mV, and parameters for the slave system are chosen with a 
difference of 10 % from the nominal values: C,„ s = 0.9 jxF/cm 2 , gK s = 32.4 mQr 1 / cm 2 , 
g Nas = 108 mQ.- l /cm 2 , g, s = 0.27 mQr x /cni-, V Ks = 10.8 mV, V N( , S = -103.5 mV and 
V/ s = —9.5517 mV. Under such a parameters both neurons shall not "share time solutions", 
then they cannot be synchronous. 
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8.3 Synchronization problem statement 



In the nonlinear control theory a synchronization problem can be stated as a feedback 
stabilization one 1611 1861 . It has been established that for a defined synchronization 
error, x e = x e ja — x e% $ (where x <G R d , e = 1,2, ...,d), there exists a synchronization 
error (dynamical) system 1611 whose trajectories exponentially converge to zero under 
a feedback control u. Hence, it can be said that the master and slave systems 
(unidirectionally) couple and attain a synchronization dynamical state under the action 
of the control command u. The following definition of Exact Synchronization is presented 
in ED: 

Definition 1. It is said that two chaotic systems are exactly synchronized if the 
synchronization error, x e = x e .M — x e ,s, exponentially converges to the origin. This implies 
that at a finite time x e .$ = x e ,M- 

In this section, a proof for the Exact Synchronization of two HH neurons is provided. 



Lemma 1. Consider the two HH neuronal systems represented in Eqs. i8.1l - X8^8v . Such 
two silent HH neurons attain dynamical states of Exact Synchronization for all t > to > 
under the action of a nonlinear feedback control despite parametric differences for any 
initial condition x,(0) = x,yif(0) — Xi t $(0) in the domain physically realizable. 



Proof. Let us define the following synchronization error system for the HH neuronal 

systems: 

x = F(x,t) + G(x)u = Af(x)+AI(t)-Bu, y = h(x)=x\, (8.9) 

where 

/ {-\/C,„ M gK M X A 2M {x XM ~VK M ) +gNa M xl M X 4M (xi,M \ 



A/(*) 



-V N a M ) +gl M ( X hM-Vl M )] + l/C ms \gK s X A 2:S {x\.s-V Ks ) 
+gNa s xl s X4,s(xi,S - V Nas ) +gl s (x LS - V, s )] } 
{0£„(xi j m)(1 - X1.m) - $n(x\ M )xiM 

-a„(xi. s )(l -x 2 ,s) + Pn{x\.s)x2.s} 

{<X,„(xi : m)(1 — *3,m) — fim(xi,M)X3,M 
-CCm(xi,s){l-X3,s)+Pm(xi,s)X3 : s} 
{ OCh (x i ,m ) ( 1 - *4,m) - Ph. (Xl )X4,M 
-(Xh(xi,s)(l -X4 :S )+ P h (x liS )x4,s} 



M(t) 



'ext M 



(t) wo 



c, 



Q 



,0,0,0 



'".V 



B = (1,0,0,0) 7 



Af(x) is a smooth vector field, AI(t ) is the difference between the external exciting forces 
and y G R represents the measured state of the system. The relative degree of a system is 
defined as the number of times one has to differentiate the output y(t ) before the control 
action u explicitly appears H58I . It can be easily shown that the relative degree for the 
system J8.9i is p = 1. Let us consider now the Proposition 4.4.2 given in |58|. The 
following stabilizing control law is proposed, 



L g Lj 1 h(x) 



^—Ljh(x) — cqIi(x) — c\Lfh(x) — ... — c p -\LPj 1 h(x) 



(8.10) 



where Lfh(x) stands for the Lie derivative of the function h(x) along the vector field /, and 
the constant parameters cq,c\ , ...,Cp_i, belong to the polynomial p(s) = cq + c\s + ... + 
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c p -\s p ~ l +s p with all its eigenvalues having negative real part. Then, the synchronization 
error system ( 18.91 is stabilized by the control law u = \—Lfh{x) — Coh(x)]/ Lq1i(x), and Cq 
is a positive real constant that represents the convergence rate. Under the control action 
given by 

U = l/C„ lM [lext M (t) ~ gK M X2M( x l,M ~Vk m ) ~ gNa M xl M X4 jM (xi,M — Vnci m ) 
~8l M ( x lM-Vl M ) } - l/C ms [lextsit) -8K s X2 >S (xi,S-Vk s ) 

~gNc, s xl s X 4 ^ (x LS - V Nas ) ~gi s (x LS - V/ s )] + C (xim - (8.1 1) 

Eq. j8.9t becomes 

x\ = -Coxi, (8.12) 

X2 = Mn(xiM)-CC n (xi t M,X\)-la n {x\ M ) + Pn(xi M )-a„(x\M,X\) 

-fin{x\,M,X\)]x2,M- [<Xn(xi,M,Xl) + P n (x LM ,Xi)}x2, (8.13) 

x 3 = Ct m (x\ M ) — a m {x\^ Ml X\) — [a m (x\^ M ) +^,h(^1,m) — Ctm{x\.M-,X\) 

-Pm(xi M ,Xi)]x3.M - [(X m {xiM ,Xi) + p,„(xi M ,Xi)]x3, (8.14) 

M = a h (x l!M )-a h (xi i M,x 1 )-[a h (xi ! u) + Ph{xi^i)-oc h (xi i u,xi) 

-f$h{xi,M,X\)}x^M- [Uh(xi,M,Xl) J rfih(x\M,Xi)]x A . (8.15) 

Because of the system (18. 1 2i -( f8~. 1 5i is minimum-phase (see Appendix A at the end of the 
chapter), the obtained control law u leads the trajectories of the synchronization error sys- 
tem to asymptotically converge to zero in a finite time. □ 

Remark 1. It should be noted that, by definition, Exact Synchronization implies that 
xm = x$ for any time t > to > and x(0) = xm (0) — xs(0) in a given domain. Now, also 
by definition, GS corresponds to the state where the states of slave systems can be written 
as a function of the master states (i.e., x$ = ^(xm))- Thus, as we shall see below, the uni- 
directional synchronization of noiseless HH neurons is exact and xs = Ixm, where / stands 
for the identity matrix. 

In order to establish the GS between the neurons ( I8.1> -( l8~8l . the following results state the 
conditions to derive an expression forx^ = \i/(xm)- Thus, we depart from the Fact 1. 
The following definitions are useful concepts presented in 1581 . For two vector fields / 
and g, both defined on an open subset U of R" (i.e., U C R"), the Lie bracket [f,g] is a 
third vector field defined by [f,g](x) = 4^/(x) — ^g(x), where ^ and -J- are Jacobian 
matrices. 

For a given set of d vector fields f\,...,fd, all defined on the same open set U, let 
A(x) = span{f\ (x), ...,/ c /(x)} be a subspace of R" spanned by the vectors fi,...,fd, at any 
fixed point x in U. The subspace A(x) of R'\ forx £ U, is called a distribution. A distribu- 
tion A is involutive if Ti G A, T2 £ A => [fi , T2] £ A, where Ti and T2 are any pair of vector 
fields belonging to A. 

Fact 1 1581 . Consider an affine nonlinear system i = f(x) + g(x)u; where x £ Q. C R", 
u £ R, g,f : R" — ► R" are smooth vector fields. Besides, let us consider that y = h(x) 
for any smooth function h(x). If involutivity condition is satisfied, then the mappings 
<t>i : R" — > R p , x h- > z and <t>2 : R" — > R"~ p , x 1— > (z,v) are such that the affine nonlinear 
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system can be written in the canonical form 

ii = Zi+i, i= l,2,...,p-l, 

z P = a(z,v) + P(z,v)u, (8.16) 
v = C(z,v), 

and can be derived from Lie derivatives of the output function h(x) along the vector fields 
f(x) and g(x) as follows 



h(x) 
Lfh{x) 



\ 



/ </>p+iM \ 

0p+2(*) 



(8.17) 



and 



(8.18) 



+1, in such a way that 
p + 1 < < n. 



(8.19) 



V = *2(*) = 

moreover, it is always possible to chose <j>f 

L g (j>j(x)=0 
□ 

The Fact 1 is well known in nonlinear control theory. Here it is included for clarity 
in presentation and exploited in neuronal synchronization towards robust feedback 
synchronization of HH neurons. 

Fact 2 lHH). If exists the map <I> = (4>i , 4> 2 ) : R" -»■ R", x i-> (z, v) derived from (EOTl and 
J8- 1 8I >. then there exists the inverse <t> _1 (<&(x)) =x € £1 C /?". This fact is proved since 
/z (jc), Lfh(x), Zy (jc) and (j> p+ i(x), (j>„(x) are linearly independent at any x in the 
neighborhood U C H C R" of the point x° in £2. 



£.4 Synchronizing the Hodgkin-Huxley neurons 



Once obtained the control action J8 .111 , it can be directly implemented in Eq. J8.5I for the 
slave system leading to the following set of coupled nonlinear differential equations, 

x i ,m = 1 / C mu [I extM (t ) - gK M x 2,M (x\. M - V Km ) 

-gNa M x\ M X AM (xi M - V NaM ) -gfa (x 1M ~Vi M )}, (8.20) 

X2,M = CC n {xi. M )(l-X 2 ,M)-Pn(xi,M)X2,M, (8.21) 

X3M = CCm (Xl Jtf) (1 - *3 M ) - Pm (xija )X3& , (8 .22) 

XA,M = &h(xijrf) (1 — Xaja) — Ph\XlJw)x4,M > (8.23) 
*1,S = l/C mM [4tf M (?) -^4,m(-«1,M- V^ m ) 

-gA^M^M-HM - ^Na M ) ~gl M {x\,M - Vfa)] 

+C (x LM -x hs ), (8.24) 

X2,S = ttn(x l . s ){\-x 2 .s)-fin{x\,s)x2,Si (8-25) 

X3.S = a m (x us )(l ~ X3 >S ) - Pm(xi,s)X3,S , (8.26) 

xas = ®h(xi,s) 0- ~ X4,s) ~ Ph(xi,s)x4,s ■ (8-27) 
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The above set of equations represents the dynamics of the HH neuronal synchronization 
when the control action u is implemented. The right-hand side of Eq. ( 18.241 describes the 
new induced dynamics of the slave system. The master and slave systems unidirectionally 
couple through u; also parametric differences have been subtracted. The term containing 
the convergence rate Co can be interpreted as a synaptic-like control current (divided by a 
constant capacitance) being CoC mM a constant synaptic conductance. In the framework of 
geometrical control and its applications on communicating systems, the synchronization 
of chaotic dynamics for the HH neurons could be understood as synchronization of a 
transmitter (master)-receiver (slave) system. Interpretation of HH neurons as chaotic 
systems in context of communicating systems, which transmit information, could provide 
insight to understand the way the biological information is processed in neuronal 
ensembles. 

Numerical simulations were carried out for the HH neuronal synchronization system. Si- 
nusoidal exciting modulations are considered, and the amplitud and frequency parameters 
are chosen within the {/-shaped curve shown in Fig. l.(b) of reference 1811 . This curve 
encapsulates the region of parameter space (in amplitude and frequency domain) where 
the exciting modulations trigger spike trains of the action potential in the model system of 
single silent HH neurons. 

Fig. 8.1 shows desynchronized regular (master system) and irregular (slave system) spiking 
patterns and the transition to a regular synchronized state of the action potentials. The 
applied forcing functions are I ex t M (t) = — 2.58sin(.245f ) and I ex t s {t) = — 3.15sin(.715f). 
Initial conditions were chosen as x,-.m(0) = (10mV,0,0, 0) andx,,s(0) = (0,0,0,0), and the 
control action was implemented at time to = 180 ms. A choice for Co = 0.3 leads to rapid 
synchronization convergence for the refractory period. The activation and inactivation 
dynamics for the ionic channels also attains synchronization state. Fig. 8.2 shows the 
evolution in time of voltage per second which is supplied to the neuronal system to 
achieve the synchronization. Fig. 8.3 shows the phase locking of the action potentials 
in synchronized state of Fig. 8.1. 




200 
Time (ms) 



Fig. 8.1: Spiking patterns of the master (solid line) and slave (dashed line) systems for 
the action potentials in desynchronized and synchronized states. The forcing functions 
amplitud and frequency parameters as specified in the text: I ex t M {t) = — 2.58sin(.245f), 

W0 = -3.15sin(.715/). 
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Fig. 8.2: Dynamical response of the implemented control action of Fig 8.1. 
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Fig. 8.3: Phase locking of the synchronized action potentials of Fig 8.1. 



Note that the nonlinear feedback (18 . 1 1 1 requires information about currents flowing 
through the membrane and the ionic channels of the master and slave neurons. Such a 
coupling cannot be implemented in practice on real neurons. However, as we shall see 
below, the controller (18. Ill allows to discuss the robust synchronization in context of GS, 
which is the most significant phenomenon in chaotic synchronization. Thus, once neuron 
synchronization is discussed in terms of GS, the robust synchronization is proposed by 
relaxing the nonlinear controller J8 . 1 11 towards a linear approach. This linear approach is 
robust in the sense that synchronization is induced in face of parameter mismatches and 
differences between amplitud and frequency parameters in external current entering into 
master and slave neurons. 
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8.5 Generalized and robust synchronization 



8.5.1 Generalized synchronization 

In this subsection, we derive the mappings <f>i : R" — ► R p , ih? and <I>2 : R" — > R n ~ p , x i— > 
(z, v), and write the HH model of neurons in canonical form (I8.16> towards generalized 
synchronization (GS). 

Dynamical models for each HH neuron can be written in nonlinear affine form i = 
f(x) +g(x)u, with 



m = 



( {l/C m [I ext -g K x 4 2 ( Xl -V K ) \ 

~gNax\xA (Xj - V Na ) - gl (Xj - V Na )\} 

a„(xi)(l -X2)-J3„(xi)x2 
a m (xi)(l -X3) -j3 m (A;i)x3 
V a ft (xi)(l-x 4 )-j3 A (xi)x 4 / 



/ 1 \ 



V J 



(8.28) 



and the output y = /z(x) is given by the membrane potential, i.e., h(x) = x\. By computing 
the Lie derivatives of the output function along the vector fields J8.28> . we obtain 



z = <& l (x) =h(x) =x\ 



and 



0i w 

= %(x) = ( 02 W 
03 (x) 

then according to Eqs. ( 18 . 1 8i and d8 . 1 91 v can be chosen as 



(8.29) 



(8.30) 







(*2 


v 2 


)- 


x 3 


V 3 1 







(8.31) 



Consequently, the previous HH model i7.4l - f7~7\ is transformed into 

Z\ = l/C m [I ext -g K v'l{zi~VK)-gNaV2V4{zi-V Na ) 

-gi (zi -Vno)] + u, 
vi = a«(zi)(l-Vi)-j3„(zi)vi, 
v 2 = am(zi)(l-v 2 )-j3 m (zi)v 2 , 
V3 = a/i (zi) (1 - v 3 ) - j3/, (zi ) v 3 . 



(8.32) 
(8.33) 
(8.34) 
(8.35) 



In what follows we show how the GS can be studied in HH neurons by departing from 
Lemma 1 and Fact 1. To this end, we can separately transform both master and slave 
neurons. In this manner, we shall derive the maps <S>m(xm) and ^(x^) to get 



z,\i 



and 



zs 
Vy 



(8.36) 



from where each HH neuron can be transformed into (I8.16> and driving signal (18. lOi 
induces the master behavior onto slave neuron. Then, if stability holds and neurons are 
minimum-phase systems, (zs,Vs) — ► (zm,v|) for t > to > 0, for any initial conditions 
(z(0), v(0)) = (4>[ (jc(0)),4> 2 (x(0))) in physical domain. Note that v| is a stable manifold 
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which can correspond to the stable manifold of the master neuron v^. In this case complete 
synchronization is achieved. In case vij ^ V^, the partial state synchronization is attained 
l59l . Anyway, the composition 4>^' (<£>] (xm); v|) = x s G i2 C /?", where i2 denotes the 
physical domain. In particular, if v| = for all time f > to > 0, where to stands for 
time of turning on the control, then xs — O^T 1 (<Pim{xm),^'2m(xm))- Since HH neurons 
J7.4ft - (l7.7i are minimum-phase systems (see Appendix A at the end of the chapter), the GS 
yields the following relation 



/ h s l (h M (x M )) \ 



xs ■ 



\ 



l 2S 

* 

K 3S 



l 4S 



(8.37) 



/ 



Now, it should be pointed out that driving signal (18. 10l > has full information about states of 
both master and slave neurons. This situation cannot be physically realizable (for example, 
currents due to the ionic channels activity cannot be available for feedback). In next 
paragraphs, a robust approach is taken from open literature to show how the HH neurons 
can be synchronized. 

8.5.2 Robust synchronization 

The nonlinear controller ( 18.1 H allows to obtain states of Exact Synchronization for the 
silent HH neurons represented by systems J8. Ii -( l8~8b . However, because of measurements 
for activation and inactivation of the ionic channels cannot be physically carried out, 
implementation of the control action ( 18. Hi would be unpractical. Then, an adaptive 
scheme to yield robust synchronization is realized by using a modified feedback control 
law. The synchronization error system ( 18.91 can be represented in the following extended 
form |6D, 



zi =ri + l3 E (z)u, f] = r(zi,Tj, v,m), v = C(zi,v), y = zi, (8.38) 

where the invertible coordinates change x = (zi, v) has been developed; rj = Afi (z\ , v) + 
A/i(f), represents an augmented state which lumps the uncertain terms (the ionic channels 
activation and inactivation variables for the master and slave systems) contained in Af\ , 
V is the state vector for the internal dynamics and Pe(z) = — 1 ■ In order to stabilize the 
synchronization error system, we consider the nonlinear controller ( 18. lit and observe that 
it can be written in the following (linearizing-like) form 



u=[r]+kzi] (8.39) 

where k E R + represents a control gain value. However, since the control law ( I8.39> 
depends on the uncertain state rj, it is not physically realizable. The problem of estimating 
(zi,?7) is solved by using a high-gain observer (dynamic compensator) 1611 . 

il = f)-«+IoKf(zi-£i), (8.40) 
f\ = L 2 kZ(zi-zi), (8.41) 

where (£i,j)) are the estimated values of {zi,T]); Lq is the unique tuning parameter, and 
represents a high-gain estimation parameter that can be interpreted as the uncertainties 
estimation rate. The parameters are chosen for the polynomial P(s) = s 1 + Ktf + K\ = 
with all its eigenvalues in the left-half complex plane. 
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The linearizing control law with uncertainty estimation that, together with the dynamic 
compensator fl8.40l -( fS~^Tt . stabilizes the synchronization error at the origin and, conse- 
quently, synchronizes the HH neuronal systems now becomes 



u=[f)+kzi}. (8.42) 



A stability analysis for the closed loop system (18.381 . J8.40i - 18.42i is provided in Appendix 
B at the end of the chapter. A tuning algorithm for stability and duration time is also 
provided in 1871 . Thus, controller J8 .40l >- J8~42l > is a general approach to synchronization 
of HH neurons despite it lacks knowledge about the states of activation and/or inactivation 
of the potassium and sodium ionic channels. 

It is pointed out that the modified feedback control law (I8.40i - d8~42l yields Complete 
Practical Synchronization 1591 . i.e., the trajectories of the synchronization error system 
converge around the origin within a ball of radius Lq 1 . 

Once obtained the modified control law, it can be implemented in systems (18 . 1 i - <!8 .81 . 
Then, we are led to the following extended system of differential equations that guarantees 
the robust synchronization of HH neurons, 

X\.M = 1 / C„ lM [I extM (t) - gK u x\jt (*l,Af ~ Vk u ) 

-gNa M 4.,M x 4M (xi,M - V NclM ) -g[ M (x lM - V[ M )] , (8.43) 

X2.M = ttn{xi.M)(l-X2,M)-Pn(xi.M)x2,M, (8.44) 

X-iM = a m (xi, M )(l-X 3 , M )-An(*l,M)-X3,M, (8.45) 

X4,M = (Xh{x 1M )(l-X4M)-P h (x LM )x4 M , (8.46) 

Xl,S = 1/C ms [l eX t s (t) -gK s X2 tS (xi,S-VK s ) 

-gNa s xl tS X4,s(xi,S-VNas) ~gl s (*1,S - V fe )] +(fl+kzi), (8.47) 

X2,s = ^(*i,s)(l-*2,s)-/i«(*i,s)*2,s, (8.48) 

X3,S = Olmfaj) (1 - X3 t s) - Pm(xi,s)x3,S , (8.49) 

M,s = (Xh(xi,s)(l-X4,s)-Ph(xi,s)x4,s, (8-50) 
zi = -kZ 1 +L Ki((x ltM -x ltS )-zi), (8.51) 
i) = I%*Z{(xi M -xif)-2i), (8.52) 



Fig. 8.4 shows the attained robust synchronization dynamics for the master (solid 
line) and slave (dashed line) systems when the modified feedback control law has been 
implemented. The control gain value was chosen as k = 1, the K\2 parameters were 
chosen for the polynomial P(s) with its eigenvalues located at s = —20, and the high- 
gain parameter is Lq — 50. The applied forcing functions are taken as in Section 8.4, 
hxt M {t) = — 2.58sin(.245f) and l e xt s if) = — 3.15sin(.715f ). Initial conditions were chosen 
as x,m(0) = (10mV,0,0, 0) and Xjs(0) = (0,0,0,0), and the modified control law was 
implemented at time fo = 200 ms. Dynamics of the ionic channels is also synchronized by 
the modified control law. Fig. 8.5 shows the evolution in time of voltage per second which 
is supplied to the neuronal system to achieve the robust synchronization. Fig. 8.6 shows 
the phase locking of the action potentials for the robust synchronization state of Fig. 8.4. 
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Fig. 8.4: Spiking patterns of the master (solid line) and slave (dashed line) systems for 
the action potentials in desynchronized state and the transition to a robust synchronization 
state when the modified feedback control law is implemented. The forcing functions are 
W(0 = -2.58sin(.245f), WO = -3.15sin(.715f ). 
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Fig. 8.5: Dynamical response of the implemented modified control law of Fig 8.4. 
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8.6 Conclusion of the chapter 



In this chapter, we have shown that there exists a nonlinear unidirectional coupling such 
that two silent HH neurons attain synchronized states in spite of parametric discrepancies. 
Results show that synchronized spiking patterns of the action potentials are displayed by 
the unidirectionally coupled system of HH neurons. The synchronization coupling via 
the control action J8.1U yields a synaptic-like control current term containing as control 
parameter the convergence rate. Increases on the control parameter have the effect to 
allow faster synchronization convergence for the refractory dynamics. Regular spiking 
patterns in synchronized state can be achieved for regular-irregular desynchronized spiking 
sequences as shown in Fig. 8.1. Because of measurements for the ionic channels activation 
and inactivation are not physically realizable, a robust adaptive scheme has been developed 
to yield the synchronization dynamics of the HH neurons. A modified feedback control 
law composed of a dynamic compensator and a linearizing control law with uncertainties 
estimation has been implemented. The adaptive scheme leads to robust synchronization 
dynamical states of the action potentials as shown in Fig. 8.4. The ionic channels activity 
is also synchronized. The dynamic compensator allows to reconstruct the dynamics of the 
states (zi,Tj) from measurements of the action potentials, and it requires only one tuning 
parameter, Lq. Artificial devices experimentally implemented in neuronal systems 1741 
are elucidating in the unveiling of mechanisms underlying synchronization and control 
parameters perform a very important role. The nonlinear control theory could provide 
useful methods in studying synchronization phenomena in single neurons and neural 
networks, even though natural properties like intrinsic noise and synaptic conductances 
must be regarded. 
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Appendix A: Internal dynamics of the synchronization error system 

The synchronization error system J8. 1 2I >-( 15TT31 is already in canonical form. This means 
that the closed-loop system has two subsystems: the first one, given by Eq. ( I8.12L 
is controllable while second one, given by Eqs. l !8.13t - i8T3l . is not affected by the 
unidirectionally synchronization force u. Hence, to assure asymptotic stability it is 
necessary to study the dynamics of subsystem J8. 13l >- J8~T5t . If subsystem J8. 13b - d8~T5l > 
is (asymptotically) stable at origin, then the closed-loop is said minimum-phase and, as 
a consequence, the synchronization force u leads the trajectories of system ( 18. 9t to zero. 
Thus, we have that, for any time t > ?o > 0, xm = x$ and GS via Exact Synchronization is 
achieved. The zero dynamics can be obtained by setting x\ = 1581 and considering the 
master system dynamics, 

x 2 = -[a„{x LM )+l5„(x 1M )}x2, 

i 3 = -[(Xm(xi&)+Pm(xiji{)]x3, 

x 4 = -[a h (x lM ) + p h (x UM )}x4, 

XI M = 1 / C mM [I extu (t ) - gK M x\ M (Xl ,M - V K „ ) 

-8Na M x l,M X4 ' M ( X l M ~ Vn "m ) ~8l M i x hM - Vl M )] ) ( A - 1 ) 

X2.M = CCn(x\,M)(\ ~ *2,m) ~ Pn(x\^l)x2M i 
X3M = CC m (xi M ) (1 —X^m) - Pm(xiM)x3M, 
X4.M = tth{x\M) (1 — X4 ! m) - Ph{xi,M)x4,M ■ 



Calculation of a linear approximation for the zero dynamics (A.l) allows to obtain the 
corresponding eigenvalues in order to determine the stability of the system. Let x z = 
f z (x z ,t) be the zero dynamics system where x z is the state vector and f z (x z ,t) is a smooth 
vector field. The linear approximation x z = Ax z where A is the Jacobian matrix of the 
mapping f z evaluated at x z = 0, yields the following result for the matrix A, 



( 



A = 



0.183 
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\ 



(A.2) 

Due to the system (A.l) has all its eigenvalues in the left-half complex plane the internal 
dynamics is locally asymptotically stable and system J8. 12l >- j8TT5l > is minimum-phase. 
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Appendix B: Stability analysis for the synchronization error system under the modified 
feedback control law 



Let e G R 2 be the estimation error vector 116 II with states given by e\ = zi — Zi and 
ei = 7] — f\. Then, the dynamics of the estimation error is represented by the following 
system, 



Because of the trajectories of the synchronization error system are contained in a chaotic 
attractor, the uncertain terms 77(f) and r(zi,TJ,V,w) are bounded functions. Moreover, 
the matrix D has all its eigenvalues in the left-half complex plane, consequently, the 
dynamics of the estimation error converges asymptotically to zero for any Lq> Lq> 0, and 
(zi j f)) — * (zi , Tj). Therefore, the closed loop system (18.381 . J8 .4-0i -( f8~42l > is asymptotically 
stable for L Q > L^> 0. 




(B.l) 
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Part III 

CONCLUSION 
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Final conclusion 



For the first part of this thesis the main original result is an efficient factorization method for 
second order ordinary differential equations (ODE) with polynomial nonlinearities. This 
method allows us to find kink particular solutions for reaction-diffusion (RD) equations and 
anharmonic oscillator equations. In addition, application of SUSYQM-type factorization 
techniques allows to find a pair of travelling wave solutions for different RD equations 
with the same wave velocity. The method is also applied to more complicated second order 
nonlinear equations with interesting results. We believe that this factorization scheme is 
easier and more efficient than other employed methods to find exact particular solutions 
of second order ODE. Exact solutions have been found for differential equations with 
applications in nonlinear physics and biology, for instance, the generalized and convective 
Fisher equations, the Duffing-van der Pol oscillator equation and the generalized Burgers- 
Huxley equation. An application to the biological dynamics of microtubules (MTs) has 
been developed as a byproduct of supersymmetric procedures. Possible interpretation of 
our results may be related to the motion of impurities along the MTs or to the structural 
discontinuities in the arrangement of tubulin molecules. Another interesting result, in the 
context of applications of supersymmetric factorization procedures in physical systems, 
is a complex parametric extension of the classical harmonic oscillator. This extension 
is based on a SUSYQM procedure that has been previously used for Dirac equation in 
relativistic particle physics. This result may have applications in dissipative (absorptive) 
processes in physical optics as well as in the physics of cavities. Also, an application to 
the chemical physics of diatomic molecules using the same supersymmetric factorization 
scheme is included. As a result an exactly solvable nonhermitic quantum Morse problem 
is obtained. 

In the second part of the thesis it was shown that two noiseless Hodgkin-Huxley (HH) 
neurons attain synchronized dynamical states when a feedback action is implemented. 
Because there exist uncertain states that cannot be accurately measured in practice (for 
instance, the ionic channels activity), a robust approach that guarantees the synchronization 
of the HH neurons is implemented. Numerical results describing the synchronized 
behavior of the membrane action potentials of the two neurons are displayed. 
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